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ABSTRACT 

A new one-dimensional, dynamical model is proposed for geometrically thin, self- 
gravitating viscous accretion discs. The vertically integrated equations are simplified 
using the slow accretion limit and the monopolc approximation with a time-dependent 
central point mass to account for self-gravity and accretion. It is shown that the system 
of partial differential equations can be reduced to a single non-linear advection diffusion 
equation which describes the time evolution of angular velocity. 

In order to solve the equation three different turbulent viscosity prescriptions are 
considered. It is shown that for these parametrizations the differential equation allows 
for similarity transformations depending only on a single non-dimensional parameter. 
A detailed analysis of the similarity solutions reveals that this parameter is the initial 
power law exponent of the angular velocity distribution at large radii. The radial 
dependence of the self-similar solutions is in most cases given by broken power laws. 
At small radii the rotation law always becomes Keplerian with respect to the current 
central point mass. In the outer regions the power law exponent of the rotation law 
deviates from the Keplerian value and approaches asymptotically the value determined 
by the initial condition. It is shown that accretion discs with flatter rotation laws at 
large radii yield higher accretion rates. 

The methods are applied to self-gravitating accretion discs in active galactic nuclei. 
Fully self-gravitating discs are found to evolve faster than nearly Keplerian discs. The 
implications on supermassive black hole formation and Quasar evolution are discussed. 

Key words: accretion, accretion discs - galaxies: active - quasars: supermassive 
black holes - methods: analytical 


1 INTRODUCTION 

Accretion discs had become a major topic of astrophysical 
research during the past five decades since the d iscovery of 
t he f i rst quasar s in_the early 1960’s dMatthews fe Sandagel 
1 19631 : ISchmidtl 1 1963h and the ir widel y accepte d theoreti¬ 
cal explanat io n by Zel’dovichl (1 19641 . fsalpeted (ll964l 'l and 
ILvnden-Behl (Il969l l. Since those days a variety of ob¬ 
jects have been identified in which mass accretion is the 
fundamental process. Among these are such different ob¬ 
jects as T-Tauri stars, catacl ysmic variables (CV) and ac¬ 
tive galactic nuclei ( AGN) llRobinsonl Il976l : iReesI If 9841 : 
lAppenzeller fe Mundtlll989l f. Thus a deep understanding of 
the accretion process forms the basis for the explanation of 
several important astrophysical processes from the forma- 
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tion of planetary systems to the evolution of super-massive 
black holes (SMBHs) in galactic centres. 


In th e _ cla ssical the o ry of disc accretion 

dWeizsackei 19481: Lust 19521 : IShakura fe Sunvaevl Il973l : 
iLvnden-Bell fe Pringlel 119741 1 viscous torques in the dif¬ 
ferentially rotating gas flow around a central object cause 
redistribution of angular momentum from the inner disc 
to outer regions. Thus the inner parts which are no longer 
fully supported by centrifugal forces move inwards in the 
gravitational potential. The same viscous stress trans¬ 
forms gravitational energy into heat leading to an intense 
radiation. 


The probably most crucial quantity regarding the ac¬ 
cretion process is the prescription of viscosity. Although 
we still lack a profound theory of the underlying processes, 
there exists broad agreement that the nature of the viscos¬ 
ity must be turbulent, at least in case of non self-gravitating 
discs. This is strongly sup ported by a simple tim e -scale ar¬ 
gument already raised by iGoldreich fe Schubertl dl967l) in 
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connection with angular momentum transport in rotating 
stars, which rules o ut any important contributi on of molec¬ 
ular viscosity (see iFrank. King fc Raind 12002 . for an ap¬ 
plication to accretion discs). The perhaps most popular 
mo del of turbulent viscosity i n this context was proposed 
by IShakura fe Sunvaevl (1 19731 ). On the basis of their vis¬ 
cosity model they were able to derive stationary solutions 
for geometrically thin accretion discs. In order to solve the 
problem they implied that the gravitational potential of 
the central object dominates thereby assuming a massless 
disc, which is in many cases a fairly good approximation. 
However, at least in AGN accretion discs there is observa¬ 
tional eviden ce that there exist dis c s with a non-Keplerian 
rotation law llGreenhill et al. 19961: ILodato fe Bert ini feoQ.'jl : 


iKondratko et al.ll2005l: |Hure et alj|201lfh~ 

Unfortunately, things become a lot more complicated 
if the mass distribution within the disc contributes consid¬ 
erably to the overall gravitational potential. In addition to 
the usual disc equations one has to solve Poisson’s equa¬ 
tion for the gravitational potential introducing a non-linear 
coupling between mass distribution and rotational velocity. 
The problem of a visco us self-gravi t ating gase ous disc has al - 
ready been tackled bv lWeizsdckeij (1 19481) an d lTrefftd (Il952l) 
who derived the basic differential equations. They discussed 
some special solutions, but did not succeed in solving the 
general problem. 

Another challenge when dealing with self-gravitating 
discs is their ability to generate stron g instabilities as was 
already pointed out bv iToomrS (11964 :). At first glance this 
seems quite desirable, because this would l ead to the pro¬ 
posed turbulence. lLaughlin fe Bodenheimeil (Il994l ) could in¬ 
deed show that simple one-dimensional diffusion models us¬ 
ing a-viscosity approximate the main properties of three- 
dimensional simulations quite well. But their analysis also 
reveals that the diffusive transport is less characterized by 
a turbulent cascade and more through the action of grav¬ 
itational torques. T his l eads t o a serious problem raised 
by iBalbus fe Papaloizoul d 19991) . They argue that a viscous 
parametrization may be inadequate because it can only de¬ 
scribe the dissipation of energy locally whereas gravitational 
forces can act over large distances and are therefore non¬ 


local. In reply to this IGammiel d200ll ) shows in a semi¬ 
nal paper that in geometrically thin discs without large 
scale structures the local treatment is applicable for sim¬ 
ple cooling models where the cooling time r c is constant. 
He derives a simple formula relating the a parame ter of the 
iShakura fe Sunvaevl viscosity prescription to f lr c . iGammiel 
furthermore shows that self-gravitating accretion discs frag¬ 
ment if t he cooling time fa lls below a critical value of 3G -1 
(see also lMeha et al.ll2005l) . 

His findings are based on the analysis of two- 
dimensional shearing bwc_simulations_aud_were_Jatercon- 
firmed by I Rice et all d2003l ) and ILodato fe Ricel d2004l ) 
in global three-dimensional SPH simulations. In case of 
more massive discs the situation seems more complex. 
ILodato fe Rice! d2005l ) report no clear evidence for global 
transport of energy induced by gravitational forces if the 

aspect ratio H[r_ _< 1/10. This was later confirmed by 

ICossins et alj d2009l ). but they show that for mass ratios 
of Mdi a c/iW* = 0.125 the fraction of non-local wave energy 
transport rises up to 1 5%. If Mdiac/M* exceeds 1/2 global 
transport dominates dForgan et al.|[201ll ). Contrary to these 


results global grid based simulations show that even in the 
low mass r egim e loca l a-models may not be applicable at 
all dMeua et al.ll200a) especially when c onsidering more re¬ 
alistic coolin g models (jBole v et al. 20061). Recent grid based 
simu lations ( Michael et al. 20121 : Steiman-Cameron et al.l 
120131) however show that averaging the results over many 
dynamic times and large spatial volumes yields ro ughly the 
radial dependence of the a parameter predicted bv IGammiel 
(l200ll) . Thus a local viscous approximation seems applicable 
as long as the disc is not too thick and too massive. 

We would like to emphasize that most of the results 
reviewed in the preceding paragraphs were obtained for 
protoplanetary discs where the aspect ratio becomes 1/10 
even for moderate disc masses. In case of AGN discs the 
situation is somewhat different, because the aspect ratios 
in these discs are thought to be smaller than those found 


in protoplanetary discs by rough l y an order of magnitud e 
(ICollin-Souffrin fe Dumonu Il990l : iLin fe Papaloizo~uril996l ) . 
Since the aspect ratio has a m ajor impact on global wave 
transport (ILodato fe Ricdl2004l ) one may model AGN discs 
with local models for mass ratios well above one. 

The modern treatment of self-gravitating accretion 
discs using sim ple one-d imensional models begins with the 
work of lPaczvnskil (Il978l ) who solves - under certain assump¬ 
tions - the vertical structure problem for geometrically thin 
discs. He also introduces the concept of self-regulation which 
is based on the idea that radiative cooling and viscous heat¬ 
ing adjust the temperature in a way that keeps the disc in a 
marginally stable state. Thereby he as sumes that turbulence 
is driv en by gravitational instabilities. ISakimoto fe Coronitil 
dl98ll) modify this work by re placing the turbulenc e mode l 
with the a-parametrization of IShakura fe Sunvaevl dl973l) . 
They derive stationary solutions and apply them to self- 
gravitating AGN accretion discs. 


Based on these early attem pts to construct self - 
gravit ating accreti o n dis c models iMineshige fc Umemural 
d 19961) and iBertinl dl997l ) found stationary self-similar so¬ 
lutions. The latter author combines the a-viscosity model 
with the concept of self-regulation (see also lBertin fc Lodatol 
Il999l) . Thereby he assumes th at the d isc is alwa ys in an 
marg inally stable state, i.e. the lToomrel parameter dToomrel 
ll964l) is close to unity. With help of this assumption he 
avoids the problem of solving the energy equation to deter¬ 
mine the temperature of the disc and thus the speed of sound 
whi ch is necessary to compute a-v iscosi ty (see S ec. l2.6D . 

IMineshige fc' Umemural d 19971 ) and iTsuribel d 19991 ) de¬ 
velop self-similar time-dependent solutions based on the a- 
prescription. The former authors additionally assume that 
the pressure scale height of the disc scales linearly with 
radius. This assumption modifies the viscosity prescrip¬ 
tion considerably because it scales with sound speed mul¬ 
tiplied by radius instead of scale height in contrast to the 
original a-viscosity. I 11 both papers the discs are isother¬ 
mal an d self-gravity is t r eated in the monopole approxi¬ 
mation. IMineshige et alj dl997l ) propose a non-isothermal 
model by introducing a polytropic relation to account for 
a fixed radial temperature gradient. They also discuss 
the possibility of Quasar formation on the basis of their 
model. All these results seem quite promising, but the time 
scales for AGN evolution always exceed the Hubble time 
for reasonable disc parameters as was a lready pointed out 
by IShlosman. Begelman fe Frankl (j 199dl ) . Another problem 
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Self-Gravitating Accretion Discs 3 


of stationary self-gravitating q-d i scs w as pointed out by 
iDuschl. Strittmatter fc Biermannl (|2000l f who demonstrate 
that one inevitably yields temperature distributions which 
are inconsistent wi t h the thin disc assumption. 

iLin fc Pringle! rtl987T ) propose a new viscosity prescrip¬ 
tion where the effective kinematic viscosity is determined 
by the typical length scale of unstable regions. They show 
that in a gravitationally unstable disc this length scale de¬ 
pends on the local mass distribution and the rotational ve¬ 
locity. With help of their new viscosity model, they derive 
time-dependent self-similar solutions. However these discs 
cannot be considered as fully self-gravitating, because the 
authors keep the central mass constant in their models and 
assume Keplerian rotation. Another problem of these solu¬ 
tions is that the viscosity model requires gravitational in¬ 
stabilities to generate a sufficiently high effective viscosity. 
Hence their effective viscosity tends to zero in the Keple- 
rian limit, because Keplerian discs are known to be stable 
(ISafronovll 195^1 . 

In a completely different approach the eff ective viscos¬ 
ity is coupled to the critical R eynolds number (IDuschl et al.l 
1 19981 : [Richard fe ZahnlFl 9991 ) . This so called /3-prescription 
relies on the observation that in laboratory experiments al¬ 
most any flow becomes turbulent at high Reynolds numbers. 
Since the Reynolds numbers in accretions discs are extraor¬ 
dinary high, one can usually expect that these flows are tur¬ 
bulent regardles s of the ac tual mechanism that generates 
the turbulence dLiist 1952lf. On the basis of this new vis¬ 
cosity prescription Duschl et al.l I 2 OO 0 I ') develop stationary 
solutions for self-gravitating accretion discs. In this context 
they also discuss the possibility of supermassive black hole 
for mation based on t heir model . 

lAbbassi, Ghanbari fe Salehil ( 2006j]_and more re- 
cently Abbassi, Nourbakhsh fe Shadmehril (l2013l l derive 
self-similar solutions for self-gravitating /3-viscous discs. 
Their models modify those fo r pol ytropic discs of 
IMineshige, Nakavama fe Umemural (Il997l ) by replacing the 
a-prescription and therefore avoid the drawbacks discussed 
above. However, by doing so they encounter a problem not 
further dealt with by the authors. In order to derive the self¬ 
similar solution they introduce a similarity variable which 
depends on the proportionality constant K and the expo¬ 
nent 7 of the polytropic relation P = Kg 1 . Both constants 
enter the set of differential equations only due to the pres¬ 
sure gradient in the radial momentum equation. As we will 
show in Sec. this term is usually negligible and the au¬ 
thors do actually neglect it by using the slow accretion limit. 
Thus, although the parameters K and 7 are removed from 
the underlying basic equations, their similarity solutions de¬ 
pend on them which causes a serious contradiction between 
model assumptions and solutions. 

Therefore, in the present paper we simplify the set of 
differential equations before applying an appropriate sim¬ 
ilarity transformation. This approach yields a single par¬ 
tial differential equation (PDE) for self-gravitating accretion 
disc dynamics (Sec. l2.5ll . Although our general derivation is 
independent of the viscosity prescription one has to select 
a specific model in order to obtain solutions of the differ¬ 
ential equation. Therefore we discuss three different viscos¬ 
ity models including the /3-viscosity ISec. 12.61) . Furthermore 
we show that for these viscosity models our disc evolution 
equation is invariant under the same scaling transformation 


which admits a similarity transformation depending on a 
single non-dimensional parameter k ISec. 13.ill . Thus we ob¬ 
tain a hole family of self-similar solutions, each with a dif¬ 
ferent value of k. We demonstrate that this parameter is 
related to the slope of the rotational velocity far from the 
origin. In addition we show that k has fundamental impact 
on the evolution of self-gravitating discs and that discs with 
an asymptotically flatter rotation law have higher accretion 
rates and therefore evolve faster than those with a nearly 
Keplerian rotation law ISec. 14.31) . 


2 THE DISC MODEL 


Our model is based on the standard theory of ge ometrically 
thin, axisym metri c accre tion discs according to IWeizsackerl 
Jl9431 and iLiistl dl952l l (for a modern treatment, see 
iKato et al.ll2008f )7 According to this we assume that the disc 
is in hydrostatic balance in the vertical direction. This allows 
us to decouple the dynamical evolution from the vertical 
structure equations by introducing the vertically integrated 
density 


E(t,r)= / g(t,r,z)dz 


( 1 ) 


and pressure 


n(t,r) = f P(t,r,z)dz. 

J-H 


( 2 ) 


Thereby the limits of integration are given by the so far 
unspecified parameter H which can be finite or infinitc0. 
We show in Sec. ED that H is related to the pressure scale 
height. It is important to mention here, that even in the case 
where H becomes large the vertical density and pressure 
gradients are rather steep, so that the thin disc assumption 
always holds. In addition to surface density and integrated 
pressure we define the vertically integrated r-y-component 
of the stress tensor: 


f H dfl 

T rv (t,r) = J t r <p(t,r,z)dz = vY,r— 


(3) 


which is usually the dominant term. Q, = v^/r is the angular 
velocity and u the kinematic viscosity, both are assumed 
to be independent of the vertical coordinate 0 in order to 
carry out the integration. The set of differential equations 
we consider in this work are then given by the continuity 
equation 



and the transport equations for radial momentum 

dv r dv r 1 <9n v 2 9$ 
~dt + Vr lfr = _ E + V ~ Ifr 

and angular momentum 


dt dt 
dt + Vr dr 


_i_d_ 

rE dr 



(4) 

(5) 

(6) 


where $ is the gravitational potential, i = rv v = r 2 Q, is the 
specific angular momentum and T rv> is given by Eq. ([jU. In 


1 For H —¥ 00 , our definition of E agrees with the usually used 
version 
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order to solve the system above one has to consider the ver¬ 
tical balance law. We derive an approximate solution of the 
vertical structure equation for an ideal gas equation of state 
assuming a polytropic relation between pressure and den¬ 
sity in Sec. 12.11 In addition one generally has to solve some 
kind of energy transport equation to determine the thermo¬ 
dynamic structure of the disc. However, we show in Sec. 12.21 
that for geometrically thin discs the system decouples from 
the energy equation if all radial gradients are moderate and 
the viscosity v does not depend on temperature. In Sec. 12.61 
we discuss some reasonable viscosity prescriptions with this 
property. 

Since the discs we examine in this work are assumed to 
be self-gravitating, we have to solve Poisson’s equation to 
obtain the gravitational acceleration ^ in Eq. 0. This is 
done using the monopole approximation (see Sec. 12.311 . 


2.1 Vertical structure 


The derivat ion of the ver tical structure equations generalizes 
the work ofjHos 3 (1977.) who studied polytropic Keplerian 
discs and |Paczvnskil ~ Jl978l 'l who also included the discs po¬ 
tential, but in a slightly different way than we do. 

The basic assumptions are hydrostatic equilibrium be¬ 
tween pressure forces and gravitational forces 


1 dP _ 
g dz dz 

and a polytropic relation according to 

P = Kg " , 0 < K, 0 < n. 


(7) 

(8) 


The constant K can be determined from the midplane values 
of density and pressure 

k = p cQ :^ = ^c? ie0e -- 

where c S)C is the midplane polytropic sound velocity 


2 

^S.c 


dP 


= ^Kq? = 2±i —. 

n ^ n „ 

=0 Sc 


(9) 


The polytropic relation allows us to express the left hand 
side of Eq. 0 in terms of g alone: 


1 dP_ _ 2 

g dz nCs,c dz\g c J dz' 


This differential equation can be integrated immediately 
using the boundary conditions g(r, z = 0 ) = Qc(r) and 
<f>(r, 2 = 0 )= $ c (r): 

/ $ _ $ \ n 

e(r ’ 2) = ec ( 1 “^rj (10) 


where the midplane values of density g c , speed of sound c s , c 
and gravitational potential 4> c are functions of radius r. The 
poly tropic index n may also depen d on r. In contrast to the 
solution given in IPaczvnsTil (j 1978l l who replaced the grav¬ 
itational acceleration in the vertical balance law 0 by an 
approximate solution of Poisson’s equation for the gravita¬ 
tional potential, our result in Eq. m is exact. However, 
since the gravitational potential <f> depends on g in self- 
gravitating discs, we cannot deduce the density profile in 
terms of analytic functions in general. 

Nevertheless, if we are dealing with thin discs, we can 


Taylor expansion around the midplane 


onal potent i; 
flLiistlll952l) 


<f>(r, z) = 3> c + z 


9$ 

dz 


z 2 9 2 $ 

~2 Ih 2 


+ O (z 3 ) . 


If we furthermore assume that the potential is symmetric 
with respect to the midplane the linear and cubic terms 
vanish and this approximation is of order 2 4 . Inserting this 
expansion into Eq. m yields an approximate expression 
of vertical density stratification for potentials with mirror 
symmetry: 


g(r,z) = - 

with the scale height 

h{r) = Cs, c 


2 n 


92£ 

dz 2 


( 11 ) 


( 12 ) 


In addition to that we define the geometric height or half¬ 
thickness of the disc as the vertical extend at which the 
density vanishes: 


H(r) = y/2 n h(r). 


(13) 


This definition of H fixes the integration limits in Eqs. 0 
and 0. In the isothermal limit n —> oo we therefore yield 
H —> oo and Eq. (Hill becomes the well known Gaussian 
3ellTi96 


profile (iLvnden-BelllllQBfll : IShakura fc Sunvaevll 1973 1 


, \ ~) 2 
£isoth {r, z) = g c e 2 V./W 


(14) 


Furthermore we may utilize the definition of H and the verti¬ 
cal density stratification gd to compute the surface density 
and with help of the polytropic relation 0 the integrated 
pressure 


V rr f~ T( n + 1) 

^T(n + |) 

(15) 

rr _ d rr {n + 1) V{n + 1) 

11 ic-tl 7T , o \ / o\ 

(n+|)r(n+|) 

(16) 


where P represents the gam ma fun c tion. These relations 
have already been derived by iHos h] dl977h in case of non 
self-gravitating discs. Dividing Eq. m by Eq. (flSll and us¬ 
ing the expression for the midplane sound velocity @ one 
obtains the useful relation 


n = rjc^c £ 


(17) 


between integrated pressure, midplane speed of sound and 
surface density. The non-dimensional function 


V 


n + 


(18) 


is always larger than 0 and becomes at most 1 in the isother¬ 
mal limit. Since rj depends on the local polytropic index n 
it will in general depend on the radial coordinate r. 

With help of Eq. (flSll and Eq. (fl3l) we can derive an¬ 
other useful equation which relates surface density to central 
density and scale height: 


E = 2\g c h. 

The non-dimensional factor 

1 7 rn T(n + 1 ) 

A — 


2 r (n + |) 


(19) 

( 20 ) 
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n 


Figure 1. The non-dimensional factor A as a function of poly¬ 
tropic index n. 


is shown in Fig. |T|as a function of polytropic index n. For 
reasonable values of n > ^ it is roughly of order one and 
becomes at most \/-7 t/2 in the isothermal limit (jLiis ti l 19521: 
iLvnden Belllll969l f . 

A nice feature of the vertical structure equations given 
above is the possibility to compute the scale height of dif¬ 
ferent contributions to the potential and combine them to 
one effective scale height: 



We demonstrate how to apply this to some important exam¬ 
ples and compare the results to those given in the literature. 


(i) point mass M*: 

GM+ 

= 


<9 2 $, 


dz 2 


GM* 




Vr 2 + z 2 ’ 

where 11 k is the Keplerian angular velocity. Hence the scale 
height of Keplerian discs is given by 


h=%£. 


( 21 ) 


This result has bee n known for decades in case of isother¬ 
mal discs (see, e.g. Prineld 198l|) an d was also obtaine d for 
poly tropic discs bv lHoshil ll 19771) an d IPaczvnskil (|l978i l . 

(ii) self-gravitating homogeneous slab: 

From Poisson’s equation we have (see, e.g. lMestelll~1963lf 


d 2 $ 

dz 2 


= A-kGqc 


z=0 


and therefore the scale height of the slab becomes 

-l 


h = c s ,cy/47rGpc 


( 22 ) 


The half-thickness computed fr om this result b y use of 
Eq. m differs from that given in lPaczvnskil d 1978T) for poly¬ 
tropic self-gravitating sheets. The deviation is of order one as 
long as n is of order one and becomes large in the isothermal 
limit n —> oo. However, in th at case our sol ution reproduces 
exactly the result derived bv ISpitzeil dl942l l. 

(iii) point mass and self-gravitating homogeneous slab: 
Combining (i) a nd (ii ) as d escribed above leads to the re- 
sult proposed bv lLris 311953) (see also lSakimoto fe Coronitl 


GM]) who derived the scale height for self-gravitating 
isothermal sheets with central point mass: 


h = c< 


c\J 4nGQc 


+ 


(23) 


Thus again our more general result is consistent with the 
isothermal limit. 

(iv) self-gravitating axisymmetric disc within an axisym- 
metric external potential $ ex t in radial balance: 

h = Cs,c(^inGg c — 2H 2 (1 + x) + A r z<I> ex t| z=0 ) 2 (24) 


where x is the logarithmic derivative of midplane angular 
velocity £2 with respect to r 


d In Q. 
9 In r ' 


(25) 


and Arz is the axisymmetric Laplacian. The derivation uses 
again Poisson’s equation for the disc potential and in ad¬ 
dition the gravitati onal balance law ( Eg. 12811 . This result 
has been derived bv lBertin fc Lodatol (jl999l l for isothermal 
discs. The contribution due to the external potential van¬ 
ishes for all r > 0 if f&ext is the point mass potential. In 
that case the non self-gravitating limit (g c —> 0, a: —> — |) 
approaches the Keplerian value of t he scale heigh t. Another 
interesting case is Mestel’s disc fsee lMestellll963l l which has 
x = — 1 and the scale height becomes that of an infinite slab 

(HU. 


2.2 The slow accretion limit 


I n this sect ion we will introduce the slow accretion limit 
(|Liist| |l952|) and show that the transport equation for ra¬ 
dial momentum (Eq. [5} simplifies to a balance law equating 
centrifugal and gravitational forces. Thereby we make use of 
some relations already derived in Sec. 12.11 The basis for the 
derivation shown below is the assumption that the rotational 
velocity of geometrically thin discs is highly supersonic and 
that the radial drift ve locity is subsonic, i.e. v r ^ c s , c v v 
(see, e.g. IPrin MEMO), in case of non self-gravitating discs 
the requirement of supersonic rotation is quite obvious and 
a direct consequence of the thin disc assumption. With help 
of Eq. (THTl one concludes 



rQ. 

Cs,c 


Vp 

Cs,c 


(26) 


We cannot derive a similar estimate if we take self-gravity 
into account, because in that case the scale height depends 
on mass distribution which in turn influences the rotational 
velocity. This makes the relation between h and v v more 
difficult. However, we will show in Section E4] that our ap¬ 
proximations are at least consistent with the assumption 
given above. 

With help of the vertical structure equation GU) we can 
eliminate n from the radial momentum equation ([5]). After 
multiplication with r/c 2 c one obtains 

r dv r f v r N 2 d In v r _ f v v \ 2 r 9<f? 

Cs ,c dt + \Cs jC J din r \c s , c ) cf )C dr 

(27) 

f | d In n cllnc^ jC d In E \ 

W n + | d In r <91nr ~*~i91nr/' 
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6 T. F. Illenseer and W. J. Duschl 


Thereby we used Eq. nil) to express the logarithmic deriva¬ 
tive of rj in terms of the logarithmic derivative of the poly¬ 
tropic index n. We can henceforth conclude that if the radial 
gradients of n, Cs,c and E are moderate, i.e. their logarith¬ 
mic derivatives are at most of order one, we can neglect 
the terms within the curly brackets multiplied by r/ ^ 1 in 
Eq. (1271) in comparison with the term (u v /c s , c ) 2 1- The 
same argument applies to the second term on the left hand 
side, which is also of order one as long as the radial drift 
velocity is subsonic. The remaining terms of the radial mo¬ 
mentum equation are 


dv r 


d$> 


r dt V<p r dr' 


In the slow accretion limit one expects that temporal 
changes of the radial drift velocity v r occur on the viscous 
time-scale r v i S 7d yn « fl _1 = r/vip. Hence we may approx¬ 
imate the term on the left hand side by 


r 


dv r 

dt 



Tdyn 

= V r V(p - <C V 

Tvis 


2 


and the radial momentum transport equation reduces to the 
gravitational balance law 


v 


2 



(28) 


This res ult has bee n know n fo r dec a des s ince the early 
works of IWeizsackeil (Il948ll and iLiis 3 (Il952l '). However, to 
our knowledge it has never been derived in such a general 
way using the vertical structure equations to rewrite the ra¬ 
dial pressure gradient. By doing so we can explicitly show 
that radial pressure forces in geometrically thin discs are 
usually rather small compared to centrifugal and gravita¬ 
tional forces. 


2.3 Monopole approximation and self-gravity 


So far we did not point out how to compute the radial grav¬ 
itational acceleration — In general there would be con¬ 
tributions from the central object and the mass distribution 
within the disc. Thus it would involve the solution of Pois¬ 
son’s equation which is difficult even in case of rotationally 
symmetric and geometrically thin systems. In this section we 
will introduce the monopole approximation for such systems 
and derive an approximate solution to Eq. (f28|) . 

_The derivation basically follows the method of iToomrel 

(1 1963h who uses Hankel transforms to compute the potential 
of razor-thin discs given in the equatorial plane by 

OO OO 

$ d (r) = -2 ttG J dkj 0 [kr) J E(s) J 0 (fcs)sds. (29) 

0 0 

Here Jo denotes the Bessel func tion of t he firs t kind of or¬ 
der 0. At about the same time iMestell (Il963l l found that 
the radial gravitational acceleration at a certain dis¬ 

tance from the centre r can be split up into three parts: 
The monopole term which depends on the enclosed mass 
(see Eq. ED and c ontributions from the mass w ithin and 
beyond r (see also iMineshige fc Umemural Il997l l. This so¬ 
lution involves spatial integrals over the mass distribution 
E(r) which cannot be evaluated analytically in general. 

However for certain centrally condensed mass distribu¬ 
tions the monopole term is dominant and the two other 


terms cancel out as was already pointed out by IMestell 
(Il963lh . Hence one may neglect all terms except for the 
monopole term and approximate the solution to Eq. (1281) 
by 


2 _ GM(r) 


(30) 


with the enclosed mass 


M(r) = M* + 2-7T 



(31) 


where M* is the mass of the central object. Since the error 
introduced by neglecting the two additional contributions 
depends on E (r) it is difficult to estimate it in general. 

In order to elucidate this we show h ow to de¬ 
rive Eq. (l30l) from Eq. (l29l) . Unlike IMestell (Il963l ~l and 
IMineshige fe Umemural dl997ll we do not split up the grav¬ 
itational acceleration obtained from (HU via differentia¬ 
tion. Instead we take the solution (12911 add the poten¬ 
tial of the central point mass and insert it into Eq. (1281) . 
Then we apply the inverse Hankel transform which provides 
us with an integral expr ession of the surface density (see 
iBinnev fe Tremainelll987l) 


E(r) = J dkJ 0 (kr) J ^v v (s) 2 — ^ kJi(ks)ds. 

0 0 

If we insert this result in Eq. ED. we can compute the 
enclosed mass. Thereby the contribution due to the central 
point mass within the integral cancels the constant term M*. 
Thus we proceed utilizing integration by parts with respect 
to s considering that 

^(Jo(fcs)) ss —kJi(ks) 

which gives us 


°° 2 

M ^ = h{f ' F ( r,s ' )d ' s ~ [%(s) 2 -^(us)] s=0 


with 


J r (r,,s)= / —Ji(kr)Jo(ks)dk = 


•F<(£) if 

•F>(5) if s>r 


■ (32) 


The solution of this definite integral depends on whether 
r < s or not. An analytical expression in terms of complete 
elliptic integrals is given in App. [Altogether with asymptotic 
expansions which allow us to evaluate the surface terms. We 
may now split the integral into an inner s < r and outer 
s > r part 


GM(r) 


= %>(0) 2 - | lim f^(s) 2 




r 


and substitute the variable of integration 


GM{r) 

r 

1 


w (°) 2 — I lim 

S—>00 


l v f( s ) 2 


1 


+ / Jr(tv(^))V<(fc)dfc-| (»*(£)) V^fcjdfc. 
0 0 
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Self-Gravitating Accretion Discs 7 


If we add 0 = +1 — 1 to the function in the first 

integral, carry out one integration and factor out the term 
v v {r) 2 we finally get the result 


GM(r) 


= v v {r) < 1 — i lim 


i u™ r f MA 


s->co S \V v (r) 


— I H(r,kr)-—+ 


f 


1 

/«(,,; 




d k 


(33) 


with 


, . v lp (s) 2 dlnv 2 o s 2 n(s) 2 / din 12 

’ u v (r) 2 dins r 2 12(r) 2 \ d In s 


+ 1 


(34) 


The comparison with Eq. © reveals that the exact solution 
of the thin disc Poisson problem reduces to the monopole 
approximation if the sum in the curly brackets of Eq. |R51) 
is close to one. This sum only depends on the rotation law 
given by 12 as a function of radial distance to the origin. 
Thus we may ask ourselves if there exist certain rotation 
laws for which this condition is fulfilled. 

First of all one should restrict the discussion to rota¬ 
tion laws for which the centrifugal acceleration r!2 2 tends to 
zero at infinity. Thus 12 oc r K with re < — | must hold in the 
limit r —> oo. Then we can neglect the surface term and the 
deviation from the monopole approximation is completely 
determined by the two definite integrals which depend on 
the rotation law through T-L multiplied by the two weight 
functions (1 — J r <)/fc and JF > /k. These functions are both 
positive and smaller than one over the whole interval of in¬ 
tegration (see Fig. IA1I in the Appendix). Hence if H does 
not change its sign, i.e. 12 is monotone, both integrals may 
at least partly cancel each other. 

The simplest monotone function one could think of as 
reasonable rotation law is a power law 12 oc r K . In this case 
one easily computes for 


/ „\ 2 («+ 1 ) 

«*(r,8) = 2(K+ !)(-) 


If we insert this into Eq. (1331) the explicit dependence on r is 
removed from the integrals. We can evaluate them numeri¬ 
cally if we specify the power law exponent re of the rotation 
law. Hence the whole term within the curly brackets de¬ 
pends only on re. This implies that M oc r 2n+3 and because 
M (r) must be a monotonically increasing function of radius 
re ^ | is required. 

In Fig. [2] the gravitational acceleration due to the 
monopole term divided by the centrifugal acceleration is 
shown as a function of re. There are two cases where the 
monopole approximation is exact: re = — | and k = —1. The 
former corresponds to Kepler ian rotation w hereas the latter 
is known as Mestel’s disc ( see lMestelll963 '). The error intro¬ 
duced by the monopole approximation would be less than 
10% for — | < re < — 1 and yet for re m — | one overesti¬ 
mates the gravitational acceleration acting inwards by less 
than a factor of 2. T his observation supports the remark in 
iLin fe Pringlel (il99Cf) that the monopole approximation is 
usually good to the 5 % level. However, as re approaches — b 
the second integral in Eq. (1331) diverges and the monopole 
approximation breaks down. 



Figure 2. Deviation from monopole approximation for power law 
rotation curve S2 oc r K as a function of the exponent. 


2.4 Supersonic rotation and self-gravity 


In Section E3 we made an important assumption for our 
model, namely that the azimuthal flow in the disc is highly 
supersonic (i>p c s , c ). We showed that this is a consequence 
of the thin disc assumption if self-gravity is negligible. Unfor¬ 
tunately it is not possible to generalize these considerations 
and simply apply them to self-gravitating discs, because v v 
couples to the surface density which is not known a priori. 

However, if the radial balance law © holds, we know 
something about the relation between matter distribution 
and rotation law. Hence we can at least check, if this relation 
is consistent with the assumption of supersonic rotation. The 
differentiation of Eq. © with respect to r implies that 


dM 

dr 


§(2* + 3) 


(35) 


where x is the logarithmic gradient of the rotation law 
fEq. 1251) . With the definition of the enclosed mass in Eq. (1311) 
we can derive another expression of its radial gradient: 


dM 

dr 


= 27t rE. 


(36) 


Hence, by equating the right hand side of Eqs. (1351) and (1361) 
we obtain: 


2ttGE = rQ 2 (2x + 3). (37) 


One can use this result together with Eq. © to eliminate 
g c from Eq. (1241) . If one furthermore substitutes point mass 
potential for external potential this finally becomes 

fe) J = i (Sl + 3 >7- 2<I + 1) (r) 1 <38) 

where A is of order unity for reasonable values of the poly¬ 
tropic index n (see Fig.Q]). This expression is the generaliza¬ 
tion of Eq. (1261) for geometrically thin self-gravitating discs 
in radial balance with gravitational acceleration given by 
the monopole approximation. It obviously simplifies to the 
non self-gravitating case in the limit x -4 —3/2. However, 
if the disc is self-gravitating, the logarithmic derivative of Q 
becomes larger than —3/2 and the term of order h/r on the 
right hand side dominates and one gets (c s , c /iv) 2 oc h/r. But 
none the less the relation (cs, c /rv) 2 <C 1 still holds and hence 
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8 T. F. Illenseer and W. J. Duschl 


gravitational balance and monopole approximation are at 
least consistent with the assumption of supersonic rotation 
in case of geometrically thin self-gravitating discs. 


2.5 The disc evolution equation 


In this section we will summarize the results obtained so 
far and derive the disc evolution equation. The basic equa¬ 
tions of thin disc evolution were given in the introductory 
paragraph of Sec. 0 In the successive subsections we showed 
that one can replace the radial momentum transport equa¬ 
tion © by Eq. J30J with the enclosed mass M(r) defined by 
Eq. 0D- Because of this it is convenient to replace the sur¬ 
face density in the continuity equation © by the enclosed 
mass as well and treat it as a time-dependent function too. 
Hence the radial integration of 0 yields 

dM 

^- = -2nrT,Vr. (39) 

With help of Eq. (1361) we may eliminate E 4 Weizsackeil 194gif : 


dM 

dt 


+ V; 


dM 

dr 


= 0 . 


(40) 


This replaces the continuity equation 0 in the set of model 
equations. 

We proceed by transforming the angular momentum 
transport equation ©. Because of Eq. (Sol we can elimi¬ 
nate v r with the quotient of temporal and radial derivatives 
of M and with help of the balance law we may replace 
M by rv^,/G in these derivatives. Thus the left hand side of 
Eq. © becomes 


dt dt _ / din M \~ 1 dt M dM 
dt Vr dr \ d In r J dt 2nT, dt ' 


In the last step we used Eq. (1361) and l = r 2 fl. If we multiply 
the angular momentum equation © by 27rrE and use the 
result above on the left hand side it transforms to 

~ rM ^ = ^^Er 2 ^) 

with the viscous stress tensor component given in Q. Again 
we can replace M by utilizing Eq. dMD and E with help of 
Eq. 437}. Thus the final result for the disc evolution equation 
becomes 

~ r4Q2 ^ = §^{ ur3Q3x ( 2x + S ))- ( 41 ) 
Equation (1411) together with the local power law exponent 
x defined in (1251) is a non-linear second order partial differ¬ 
ential equation which describes the advection and diffusion 
of angular velocity under the influence of self-gravity and 
viscous friction. 

If the kinematic viscosity i/(r, f) is given as a function 
of radial distance and time, one may in principle derive a 
solution to the disc equation provided that one specifies ap¬ 
propriate initial and boundary conditions. Even if v is not 
given explicitly, the equation remains solvable if the viscos¬ 
ity depends on directly or indirectly through E or M. We 
will discuss some possibilities for an appropriate viscosity 
prescription applicable to accretion discs in the next sec¬ 
tion. _ _ 

The disc evolution equation derived by iTrefft d l| 19521 ) 
seems to be quite similar compared to our equation. How¬ 
ever, we would like to emphasize that there exists a very 


impor tant difference. In contrast to our approach ITrefft j 
4195211 solves Poisson’s equation for the disc potential in a 
pure two-dimensional world. Therefore she obtains Mkd? 
as the radial balance law. This result differs considerably 
from the solution we yield for an infinitesimally thin disc 
embedded in a three-dimensional space. 


2.6 Viscosity prescription 

A proper description of the viscosity coefficient v is a peren¬ 
nial problem when modelling accretion discs. Although there 
has been a long lasti ng debate since the early works of 
IWeizsackeil 4l948l l and lLiistl 4l952l l on the nature of the vis¬ 
cosity coefficient it became a broad consensus that molec¬ 
ular viscosity is not su fficient to explain the accretion pro¬ 
cess 4Frank et alj l2002f ). In case of non self-gravitating discs 
there are basically two processes that lead to sufficiently high 
stresses: Strong, large scal e magnetic fields and turbulence 
JShakura fc Sunvaev Il973l f . 

However, in self-gravitating discs there exists a com¬ 
pletely different mechanism for redistribution of angular mo¬ 
mentum. If these discs become gravitationally unstable, they 
can transfer angular momentum via compressible density 
wave s and dissipate energy in shocks iBalbus fc PapaloizorJ 
Il99fllf . Since the amount of angular momentum transfer de¬ 
pends on the details of the inhomogeneous structures, one 
usually has to perform multidimensional simulations to in¬ 
vestigate these processes. Thus it has been the fundamental 
question of self-gravitating accretion disc theory in the past 
20 years whether or not it is possible to model those discs 
with a simple one-dimensional diffusion model. Since we do 
not want to repeat ourselves we refer to Sec. |T] for a short 
review of the most important findings including important 
references. We proceed summarizing the discussion on this 
issue with the statement that the viscous diffusion approxi¬ 
mation seems possible even for self-gravitating discs if they 
are geometrically thin and not too massive. 

In the context of non self-gravitating discs a very pop- 
ul ar description of the effectiv e shear stresses was proposed 
by IShakura fc Sunvaevl (l973j) who derived a parametriza- 
tion which couples the dominant stress tensor component 
tr V to pressure 

tnp = —aP (42) 


with model parameter 0 < a < 1. If we carry out vertical 
integration according to Eqs. © and ©, we can cast this 
equation with help of Eq. 113 i nto an expression for the 
effective turbulent viscosity fsee lLin fc Pringl3lll990l 7 PI 

< 43 > 

The new parameter a = —arj/x is again smaller than one, 
because rj ^ 1 (see Eq. (1181) 1 and the log arithmic derivative 
of il ( Eq. 1251) is negative and of order one. IBalbus Hawlevl 
4l991i) have discussed an instability that in magnetic accre¬ 
tion discs can give rise to viscous stresses of the a type with 
if only marginally - the required strength. 


2 In case of negligible self-gravity one can utilize Eq. (1261) to 
transform the viscosity prescription into the well-known formula 
v — ahc s ,c- 
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Self-Gravitating Accretion Discs 9 


The a viscosity couples to the midplane speed of sound 
c s ,c which in turn depends on midplane temperature. There¬ 
fore it is generally not possible to compute the viscosity coef¬ 
ficient without knowing anything about the temperature dis¬ 
tribution within the disc. Unfortunately this would require 
solving the energy equation in addition to the disc ev o lution 
equation derived in the previous section. |Pac.zvnsl3 dl978l 'l 
was the first who came up with the idea of self-regulation 
(see also lBertinlll997l : lBertin fc Lodatolll999l : lLodatoll2007h 
which circumvents the solution of the energy equation by 
simply proposing that the flow in self-gravitating discs is al¬ 
ways at the b order of instab ility Q « 1 with the Toomre 
parameter (see lToomrelfl964l') 


Cs,cAe ^ 

7rGE ~ ttGE 


(44) 


where K e is the epicyclic frequency 


K e = Q, 


<9 In l 2 
i91nr 


Qy/2x + 4 


fl. 


The error induced by the approximation is of order unity 
for any reasonable value of *. If we solve Eq. (l44ll for E and 
insert the result into Eq. (EH) we obtain an expression for 
c s ,c in terms of v v and Q: 


Cs,c = (x + |) Qv v 


(45) 


which allows us to eliminate Cs, c from the viscosity prescrip¬ 
tion (1431) . Hence we finally get 

v = j3(2x + 3) 2 r 2 0 (46) 

with the new parameter /3 = aQ 2 /4. We would like to men¬ 
tion that the assumption of marginal stability, i.e. Q ~ 1, 
might be a problem, because this would contradict the 
assumption of supersonic motion for self-gravitating discs 
where the logarithmic derivative of the rotation law deviates 
from its Keplerian value x = —3/2 (see Eq. (14511 and discus¬ 
sion in Sec. 12.411 . The viscosity p r escrip tion in Eq. (1461) was 
first proposed by iLin fe Prin glj dl987l . hereafter LP) who 
assumed that the turbulent viscosity is caused by some kind 
of gravitational instability: 


where L cr i t is the typical length scale of unstable modes. If 
we eliminate E using Eq. S3 we can transform this equation 
into 


/ 2® + 3 
V 2tt 


r 2 n 


which becomes exactly (1461) if we set /I = 1/ 47T 2 ~ 0.025. It 
was alre ady mentioned bvlBertin fc Lodatol dl999l f and oth¬ 
ers isee ILin fc P ring le 1 19901) that LP viscosity prescription 
fails inevitably in the Keplerian limit, because self-regulation 
can only occur in the self-gravitating regime where the disc 
can become gravitationally unstable. The consequence is 
that u vanishes as x —> —3/2. The reason for this odd be¬ 
haviour is the assumption that one can assign a finite and 
fixed value to Q and therefore /3. However, in the Keple¬ 
rian limit the disc must be Toomre stable and therefore Q 
must become much larger the one. This could compensate 
the factor (2a;+ 3) 2 in the viscosity prescription which tends 
to zero as x —>■ —3/2. To shed more light on this, we insert 


the ratio c s , c /iV from Eq. (1451) into Eq. (1381) . If we assume 
that the ratio h/r is small, but always larger than zero in 
the Keplerian limit as x —> —3/2, we conclude that the left 
hand side of 

(x + §) 2 Q 2 = (2x + 3) ^ — 2(a; + 1) 


must remain greater than zero too and therefore the vis¬ 
cosity should not vanish if self-gravity becomes negligible. 
However it might become rather small, because h/r <C 1 in 
thin discs. 

The simplest modification of (1461) to overcome these lim¬ 
itations would be to ass ume that v oc r 2 Q. A visco sity of this 
kind was proposed by iDuschl et alj (Il998l . l200d . hereafter 
DSB) who assumed that the effective Reynolds number does 
not fall below the critical Reynolds number. The so-called 
/3-ansatz just reads 


v = /3 rvtp = /3 r 2 fi. 


(47) 


where the constant parameter (3 is given by the inverse of the 
critical Reynolds number. This is still a parameterization 
and not an ab-initio solution to the problem. The ansatz 
S3, however, has several promising aspects: 


(i) For fully self-gravitating disk regions, which are dom¬ 
inated by the mass distribution of the disk in the direc¬ 
tion perpendicular to the rotational plane as well as in the 
rotational plane itself, ansatz (S3 is equivalent to the hy¬ 
pothesis that the ratio of dynamic timescale (rd yn ~ fl -1 ) 
and viscous timescale (r v i s ~ r 2 /v) is constant. In other 
words one assumes a linear relat ion between the two rel - 
evant timescales: r v i B = /3 _1 ra yn . IDuschl fc Britschl (120061 1 
succeed in showing numerically that instabilities in a self- 
gravitating flow lead to a viscosity of the functional form 
of the f3- ansatz, indeed. This is furthermore supported by 
the observation that in the fully self-gravitating limit, where 
the rotation law deviates considerably from Keplerian mo¬ 
tion (a: > —3/2), the /3-prescription approaches the values 
obtained with the self-regulated a-ansatz. (S3). 

(ii) For Keperian self-gravitating disk regions, the vertical 
structure is dominated by local self-gravity, i.e. , the local 
mass distribution in the disk, while in the rotational plane 
the (almost) equilibrium between gravitational and centrifu¬ 
gal forces is still determined by the central mass, i.e. , the 
rotational velocity is Keplerian. There a functional form of 
S3 solves the p roblem of a quasi-thermostat as discussed 
by Duschl et all d2000lf . 

(iii) In the limit of negligible self-gravity (Keplerian 
disk regions) it smoothly merges into the a-prescription 
(IDuschl et alJ Il998l : IDuschl et al.l l2000l l . This can be seen 
by plugging the non-self-gravitating (Keplerian) scale-height 
m in the /3-viscosity prescription S3: 


/3— = /3 f-V 

p n p \h) n 


A comparison with the a-ansatz (Eq. 1431) reveals the above 
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10 T. F. Illenseer and W. J. Duschl 


mentioned relation^. Here, /3 is equivalent to a constant 
Reynolds number. 


In that respect the /3-ansatz is, at least formally, a general¬ 
ization of the classical a-ansatz, which otherwise runs into 
trouble in the fully and Keplerian self-gravitating domains. 
To couple all three regions, /3 needs to be chosen such that a 
smooth transition between them can be achieved. If, similar 
to the reasoning for classical accretion disks, one takes re¬ 
sort to time scale arguments, it turns out that the standard 
values of a, which are of the order of 10~ 2 " are compatible 
with values of /3 which correspond to those of the inverse of 
the critical Reynolds nu mbeQ Some asp e cts of this have 
also been investigated bv lRichard fe Zahnl (Il999l . hereafter 
RZ) who proposed a quite similar prescription: 


v = /3 r 3 


dQ 

dr 


= /3 |rc| r 2 fi. 


(48) 


It differs from the previous one by a factor of |x| which is of 
order unity in most cases. 

All the parametrizations described above differ only in 
their dependence on the local power law exponent x. Hence 
we may combine the three viscosity functions into one for¬ 
mula and introduce a new function f{x ) to distinguish be¬ 
tween these prescriptions 


f 1 DSB 

v = /3r 2 flf(x) with f(x) = < |x| RZ (49) 

[ (2x + 3) 2 LP 

The magnitude of the viscous coupling constant (3 proposed 
by LP and DSB is of order 10 -2 to 10 -3 whereas RZ derived 
smaller values of approximately 4- IIP 6 . However, the precise 
value of /3 is not of substantial importance for the solution 
of the disc evolution equation. If we define the new time 
variable r = (3 1 we can rewrite Eq. (S3 thereby eliminating 
P from the equation 

- = ^(r 5 n 4 x(2i + 3)/(a;)). (50) 

Hence we can transform any solution f i(r, t) back to fl(r, t) 
by rescaling the time variable. Although Eq. (1501) does not 
depend on /3, its actual value may have an impact on the 
solution due to boundary conditions (see Sec. 13.3.21) . 


3 SELF-SIMILAR SOLUTIONS 

In this section we show how to solve Eq. (ToUl) using sim- 
ilarity method^. To simplify the derivation we introduce 

3 The aspect ratio h/r in Keplerian discs usually depends only 
weakly on r. It is roughly of the order of 10“ 1 . .. 10 —2 in pro¬ 
toplanetary discs llAndrew s et al.| |2009|) and an order of mag¬ 
nitude smaller in AGN discs ICohin-Souffrin A Dumond 1 1 99(1 : 
iLin fc Paoaloizoull 1 99fi ). 

4 Despite the fact that in a gravitationally driven disk, the 
Reynolds number loses its meaning, a number with a similar func¬ 
tional dependence shows up in the self-gravitating context, too, 
as the ratio of viscous and dynamical time scales, r v i s and tj yn , 
respectively. In that respect, assuming a constant value of f3 is 
tantamount to a constant ratio of the two time scales. 

5 For an elaborate discussion of Lie groups, similarity meth- 
ods and related t opics the reader may consult the textbook by 
iBluman &; A ned d2002l ). 


non-dimensional variables and functions and rewrite the disc 
evolution equation in terms of these. If we specify an arbi¬ 
trary length scale r and mass scale M and define the time 
scale f according to 



we can non-dimensionalize all variables, functions and equa¬ 
tions derived in the previous sections. The equations re¬ 
tain their form except for those containing Newton’s grav¬ 
itational constant G which must be set to unity in the 
non-dimensional equations. Once we have solved the self- 
gravitating accretion disc problem for the non-dimensional 
functions we can use these scaling parameters to switch 
back to the real world quantities. In Sec. 14.51 we discuss the 
self-similar evolution of massive AGN accretion discs and 
demonstrate how to apply the inverse transformations to 
recover the dimensioned quantities. 


3.1 Differential equation of self-similar evolution 

A short calculation yields that Eq. dEO) is invariant under 
the family of one-parameter Lie groups of scaling transfor¬ 
mations 

r = A “r, r = A V, $Y = A (52) 

if and only if c = —b. The group invariant^)] are 

£ = rr 1,/K and y = — (kOt) -1 (53) 


with k = —b/a the general group invariant solution is de¬ 
termined by the expression F(y,£) = 0 where F is a (not 
yet determined) function of the group invariants. This is an 
implicit definition of the function y(£). Thus with Eq. (1531) 
one can write down the explicit form of the group invariant 
solutions 


0.(r,r) 


1 

Kry(^(r,r)) ' 


(54) 


The remaining problem is to determine y as a function of 
the similarity variable £. 

If we apply the original PDE (1501) with x given by 
Eq. (1251) to the group invariant solutions we obtain a cou¬ 
pled system of ordinary differential equations (ODEs) for 
*(£) and y(£) 


dx 
d In £ 
d y 

d In £ 


(x — K)y — (4a: + 5 )z(x) 
z'(x) 

= -xy 


(55) 

(56) 


where z (x) is the derivative with respect to x of the viscosity 
dependent function 


z(x) = x(2x + 3 )f(x) 


( x(2x + 3) DSB 

< —a: 2 (2a; + 3) RZ (57) 

[a:(2a: + 3) 3 LP 


and f(x) is given by Eq. (1491) . In case of the RZ-prescription 
we assumed that x < 0, i.e. Q always decreases as r oo. 


6 These assignments are ambiguous because any function of the 
group invariants is again a group invariant. 
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The remaining constant k in Eq. (1551) is a free parameter 
which has to be determined from the auxiliary conditions 
(see Section m - 

The planar differential system (1551561) is an autonomous 
set of two first order ODEs. It is always possible to reduce 
these systems by elimination of the independent variable to 
a single first order equation 

dx (4x + 5)g(j) - (x - k)y . . 

d y xz'(x)y 

The solution of Eq. (1581) gives x(y) which can be used to 
integrate Eq. (l56l) 


hi£ 



(59) 


This yields £(v) an d its inverse y(£) which is required to 
compute the self-similar solution f2(r, r) given in Eq. (1541) . 
Unfortunately, Eq. (1581) is a non-linear ODE and one can¬ 
not expect to solve it in terms of known functions, in gen¬ 
eral. In fact, the inverse of Eq. < 1551 ) is an Abel Equation 
of the second kind. We were not able to identify a non¬ 
linear transformation which maps our equation to any of 
the few classes with kno wn analytical solutions (see, e.g. 
IPolvanin fc Zaitsevll2003l '). 

Although analytical solutions are not excluded per se, 
we focused on numerical methods to solve the self-similar 
disc equations. The task is basically to solve the autonomous 
system of ODEs (1551561) which can be achieved using stan¬ 
dard techniques. However, the system under investigation 
exhibits some pitfalls such as singular points where the so¬ 
lution is not unique. Careful treatment of these peculiarities 
is essential. We discuss the numerical solution procedure in 
detail in Sec. m 


3.2 Auxiliary conditions 

The problem addressed in the previous sections is a so-called 
initial-boundary-value problem. This problem is well defined 
only if we provide auxiliary conditions in addition to the 
PDE (1501) describing the dynamical evolution of the func¬ 
tion Sl(r, r). Thus we have to define some initial state flo(r) 
at time r = 0 from which the solution evolves and we must 
specify how Q changes on the boundary of the spatial do¬ 
main for all r > 0. Since the spatial domain we are interested 
in is the positive real axis, its boundary is determined by the 
limits r = 0 and r —> oo. Therefore one has to specify three 
auxiliary conditions. 

If we are looking for self-similar solutions, we demand an 
additional restriction to the set of possible solutions which 
must be related to the auxiliary conditions. The self-similar 
solutions are obtained from the integration of the first or¬ 
der system of ODEs ED 1561) . In contrast to the original 
PDE, this problem is well posed if we supply initial condi¬ 
tions for x and y. Hence the number of necessary auxiliary 
conditions is reduced by one. As a consequence self-similar 
solu tions only ex ist if two of the original conditions coalesce 
fsee lAmeslll965l . Chap. 4.3). Thus, although one can nor¬ 
mally choose arbitrary auxiliary conditions, the requirement 
of self-similarity imposes some restrictions (see Sec. 13.2.21) . 

In addition, one has to ensure that these conditions are 
consistent with the underlying physical model. In case of 
the disc evolution model described above there are some 


constraints due to the fact that the surface density E must 
be positive everywhere at any time. The same applies to the 
enclosed mass which is not only positive but also a mono- 
tonically increasing function of radial distance d r M ^ 0 for 
all r > 0 (see Eq. (1361) 1. Since D is related to M via Eq. (1301) 
the auxiliary conditions (and the solutions too) must fulfill 

r 3 fl 2 >0 44 D ^ 0 (60) 


and 


dInr 3 P 2 

cllnr 


= 2x + 3 > 0 


44 x > -§ 


( 61 ) 


for all 0 < r < oo and 0 ^ r < oo. Condition (16011 may 
become even more restrictive if we demand that D must be 
continuous. Therefore we conclude that 


D > 0. 


(62) 


The sign is determined by the orientation of the rotational 
axis, which we define to point always in the positive 2 - 
direction. 


3.2.1 Initial conditions 

Let us suppose that the initial condition flo(r) satisfies the 
relations (I6T1I621) . With Eqs. (l53ll we obtain 

L2o(r) = - lim (KTi/(£(r,T)))~ 1 = -KrMim 

T—> U £—>00 

where we assumed that k < 0 (which will be justified below). 
The limit on the right hand side is independent of r and 
should be finite. Hence we conclude that 

y(£) oc as £ —> oo (63) 

which has the implication that the initial condition for self¬ 
similar solutions must be a power law of radius 

Do(r) oc r K . (64) 

In view of this result the assumption that k < 0 seems rea¬ 
sonable because otherwise the initial Do would be an in¬ 
creasing function of radius which we want to exclude from 
our considerations (see Sec. IP1) . We therefore conclude that 
the so far unknown parameter k in Eq. E3 introduced by 
the requirement of group invariance imposed on the solution 
is simply the power law exponent of the initial condition. 

We already discussed in Section [2.31 that rotation laws 
D oc r K with an exponent k greater than — I cause infinite 
centrifugal forces as r —» oo. Furthermore we showed that 
the monopole approximation breaks down if the power law 
exponent approaches — \. Taking into account that Eq. ED 
holds, one should therefore demand that 

-§<«<-!• (65) 

This restricts the parameter k to a very limited range of ac¬ 
cessible values between centrally condensed mass distribu¬ 
tions and configurations with flattened rotation curve where 
the mass is initially dispersed over the whole disc. 


3.2.2 Boundary conditions 

Normally one can think of a variety of physically meaningful 
boundary conditions for the self-gravitating accretion disc 
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problem. However, in case of self-similar solutions the se¬ 
lection of valid boundary conditions is restricted, as was 
already mentioned above. Since £ depends on r and r ac¬ 
cording to (53J we have 

{ r —> 0 for any fixed 0 < r < oo 

r —> oo for any fixed 0 < r < oo 

if k < 0. Therefore the outer boundary condition must coa¬ 
lesce with the initial condition. 

This result becomes quite clear, if one remembers that 
the dynamic time scale for the discs evolution is roughly 
given by the inverse of St. Since the initial condition is a 
decreasing power law, its inverse tends to infinity as r ap¬ 
proaches infinity. The evolutionary time scale becomes in¬ 
finitely large which means that the disc does not evolve at 
all and stays in its initial state at the outer rim. 

At the inner boundary there are basically two types of 
reasonable boundary conditions: Vanishing and finite torque 
supplied at the origin. The viscous torque G(r, r) is given 
by 

G(r, t) = 2nr 2 T rv (66) 


3.2.3 Conservation conditions 


Instead of the initial condition one could also specify the 
required auxiliary condition using a different but equivalent 
formulation. Let’s define 

rR(r) 

C(t) = 2n Efi'r rdr 

Jo 

with arbitrary constant exponents a and 6. The quantity C 
is conserved if 



(69) 


For example if a = 0 and 6 = 0 then C is the enclosed mass 
within radial distance 


R(r) = £o t 1/1 


where £o is a fixed value of the similarity variable £ fEg. 1531) . 
We may rewrite the definition of C using non-dimensional 
variables and functions: 


fo° (2x + 3)y- (a+2) t b+2 dt 

^•(a+2)^-(a+2)-|-(b+3)/rt; 


(70) 


where the stress tensor component T rlp is defined in Eq. ©. 
One easily verifies that this could be rewritten in terms of 
*(£)> J/(£) an( i the similarity variable £: 

G(r, t) = /3fiT 4 T _ ( 5 / K + 4 ) £ s y~ 4 z(x). (67) 

For any fixed and finite time r Eq. (1531) yields 


Since C must be finite in order to be conserved, the integral 
in the nominator should give a finite number and condition 
(Eg) is fulfilled if the exponent of r in the denominator van¬ 
ishes, i.e. 


6 + 3 
a + 2 


(71) 


We therefore conclude that the torque vanishes at the inner 
boundary at any time 0 < r < oo if 

lim £ s y~ 4 z(x ) = 0. (68) 

«->o 

An alternative would be that this limit is finite. In that case 
there are three distinct ways how the torque acts on the 
inner boundary depending on the value of «, namely 

k < — | the torque decreases with time 

k = — | the torque remains constant 

k > — | the torque increases with time. 

Therefore a constant torque boundary condition is only ap¬ 
plicable if the initial condition has a unique slope. We would 
like to emphasize that this restriction is due to the require¬ 
ment of self-similarity imposed on the solutions and not a 
limitation of our disc model. However, one should not be 
worried too much, because in any realistic accretion disc 
scenario one would not expect that the torque acting at 
the inner boundary remains constant over a significant time 
span. In fact the most reasonable assumption in case of discs 
around black holes seems to be the zero torque boundary 
condition, because the event horizon prohibits any coupling 
to the central object. Another possibility would be a central 
object of finite size like an inner - possibly geometrically 
thick - accretion disc or a spinning protostar with spatial 
extent much smaller than the typical scales of the surround¬ 
ing self-gravitating disc. In that case the torque acting at 
the inner boundary may indeed be finite and it is most likely 
that it decreases with time as the central object spins down. 


The initial power law exponent r is completely determined 
by the numbers a and 6 which define the conserved quan¬ 
tity C. Therefore one concludes that each initial condition 
generates a self-similar solution with specific conservation 
properties. 

For example, if a = 0 and 6 = 0, C is the enclosed disc 
mass within radius R(t) which is conserved if k = — |. This 
is a reasonable result because in that case the mass of the 
disc must be finite and therefore the rotation law at the 
outer boundary as r —> oo is necessarily Keplerian. We will 
discuss these solutions in more detail in Sec. EH 

Another interesting case is a = 1 and 6 = 2 where C 
corresponds to the discs angular momentum. Unfortunately, 
this leads to k = — | which conflicts with the requirement in 
Eq. J65J. Thus a physically meaningful self-similar solution 
with conserved angular momentum within the disc through¬ 
out the whole evolution does not exist. This does not mean 
that all other solutions violate angular momentum conserva¬ 
tion. Instead, these discs do exactly what they are supposed 
to do, namely they transfer angular momentum from the 
inner regions to the environment beyond the radius R(t). 


3.3 Numerical method 


As already mentioned in Section l3TTl the initial value problem 
defined by the system of non-linear ODEs (Eqs. [55] and 1561) 
can be solved numer ically with standard techniques . We 
use the p rograms o de ( Tufillaro, Abbott fc Reillvll 19921 ) and 


lsode rtHindmarshlll983h which implement different numer¬ 
ical integration schemes including explicit and semi-implicit 
methods such as Runge-Kutta-Fehlberg and predictor- 
corrector Adams-Moulton schemes as well as fully implicit 
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Figure 3. Phase diagram of the planar differential system defined 
by Eqs. (15515611 for DSB viscosity prescription with k, = —6/5. 
The arrows show the local gradient of the integral curves with 
direction pointing towards increasing values of the parameter £ 
on the curve. Some selected solutions obtained for different initial 
values are shown as solid curves. 


schemes like Gear’s method wh ich utilizes Backward Dif¬ 
feren tiation Formulas (BDF) (see lHairer, N0rsett fc Wanned 
1 19931. and references therein). 

For comparison we tried another numerical scheme 
which solves differential algebraic equations (DAE). Thereby 
we must transform the system of ODEs into a system of 
DAEs which was simply done by adding z(£) to the set 
of independent functions. Then rewriting Eq. ESI) yields 
a differential equation for z(£) which must be combined 
with Eq. ESJ) and the algebraic constraint for a;(£) given by 
Eq. (1571) . As for ODEs there is a variety of different numer¬ 
ical methods available to solv e DAEs ( see lHairer fe Wanner! 
1 19961 : iPetzold fc Asche rlll998lf . We used the prog ram daspk 
described in Brown. Hindmarsh fc Petzoldl d 1994 1. 

We found that the results obtained with the different 
programs and numerical schemes are similar within limits 
of numerical errors. In view of performance and numerical 
efficiency we would certainly recommend the BDF scheme 
and the DAE solver which allow for larger step sizes. This 
is important for large £ because the system becomes stiff in 
this limit. 

Moreover, in order to solve the problem one has to de¬ 
termine feasible initial values for the numerical integration. 
This basically implies that the solutions obtained with these 
initial values must satisfy the restrictions and the auxiliary 
conditions discussed in the previous section. But there are 
still some free parameters, namely the proportionality con¬ 
stant in the initial condition fEq. 1611) and the torque applied 
to the disc at the inner boundary (Eq. 1671) . 


3.3.1 Phase Plane 

In this section we examine the phase plane of the planar dif¬ 
ferential system to get an idea of the solutions topology. By 
doing so we can check which part of the parameter space is 
admitted by the restrictions given above. Since the differen¬ 
tial equations depend on the parameter k and the viscosity 


Figure 4. Phase diagram around the improper node for DSB 
viscosity prescription with re = — 1; values on the y-a.xis are scaled 
by a factor of 1000. 


prescription through the function z(x) (Eq. 1571 1 the phase 
diagram shown in Fig. [3] would change if these were varied. 
However, the main features remain unchanged. Among these 
are 


• two nodes on the z-axis at x =$ — | and x = 0 

• one saddle point on the aj-axis at x = — (j 

• one singular vertical line. 


The location of the singular line x c depends on the viscosity 
prescription (—| for DSB, —1 for RZ and — § for LP). In 
addition, there exists another singular point - except for the 
case where k = x c - which is always located on the singu¬ 
lar line. On the singular line the numerical integration fails, 
because the denominator on the right hand side of Eq. ESI 
vanishes which causes the derivative of x with respect to 
ln£ to become infinitely large as x —> x c . The only way 
a solution may continuously pass it, would be through the 
fourth singular point where the numerator vanishes simulta¬ 
neously. The type of this singular point and its location on 
the singular line depends on viscosity prescription and re: 


DSB: 


9 

4k + 3’ 


RZ: 


1 

K + 1 ' 


T r) 15309/128 
1 8 k+ 3 


The point disappears (its y-value tends to infinity) in all 
cases if k equals the location of the singular line x c . It lies 
below the z-axis if k < x c and otherwise above it. 

In Fig. El the shaded area marks the region where the 
solutions match the requirements ED and ED- In addition 
we set an upper limit on x to exclude rather flat rotation 
curves where the monopole approximation no longer holds 
(see Sec. 12.31) . There are basically two distinct families of 
integral curves of interest. One starts somewhere on the sin¬ 
gular line x c > k and approaches y —> +oo as x —> k + 
from above and the other originates at the singular point 
(—§,0) and approaches y —> +oo as x — > re - . Both fami¬ 
lies exhibit similar asymptotic behavior for large y which is 
compliant with the required initial and outer boundary con¬ 
ditions, but only the second family admits solutions with 
the correct asymptotic behavior as y —> 0 namely x —> — 4, 
i.e. Keplerian motion as r —> 0. 


© 0000 RAS, MNRAS 000 , 000-000 




















































14 T. F. Illenseer and W. J. Duschl 


Further analysis therefore focuses on the shape of the in¬ 
tegral curves in the vicinity of the singular point at (— |, 0) 
which corresponds to the inner boundary of the accretion 
disc. In case of DSB and RZ viscosity this singular point 
is hyperbolic and can be classified as an improper node 
(see Fig. 0- Its characteristic directions are the x-axis along 
which infinitely many solutions approach the singular point 
and the straight line 


»(*) = (I) 3 ' 



with 



for DSB 
for RZ. 


(72) 


Along this line there exist exactly two distinct solutions ap¬ 
proaching the singular point from the two opposing direc¬ 
tions. If k = — | it becomes the vertical line x = — |. 

The characteristic directions separate the reasonable so¬ 
lutions from those which contradict the requirements men¬ 
tioned previously. Therefore only solutions from the upper 
right quadrant which pass a point between the characteris¬ 
tic directions (shaded region in Fig. 0 will proceed to the 
singular point with y > 0 and x > — §. The analysis of the 
linearized problem reveals that 

y(x) oc (x+ |) 5 


for all solutions approaching the singular point along the 
x-axis whereas the solution along the other characteristic 
direction follows the straight line given by Eq. m- All these 
solutions obey x —> — | as y —> 0. Since —x is the local 
power law exponent of ?/(£) (see Eg. 1561) one concludes that 

3 . 2 

y oc £2 for small y or inversely £ oc as £ —0. Hence we 
can compute the limit in the boundary condition given by 
Eq. m 

lim £ 5 y ~ 4 z(x ) oc lim (x + |) 9 z(x) 


where q = — | for the distinct solution which follows the 
straight line and q = — 1 for all other solutions approach¬ 
ing the singular point along the a:-axis. We infer from the 
definition of z(x) in Eq. (15TI) that the limit is zero (vanish¬ 
ing torque) only for q > — 1 and finite for q = — 1. Hence 
there exists a unique solution which matches the no-torque 
condition at the inner boundary, namely the solution with 
q = — Moreover, there are infinitely many solutions each 
of them associated with a distinct function describing the 
time dependence of the torque at the inner boundary. 

If we change over to the LP viscosity prescription we 
observe that the shape of the integral curves near the singu¬ 
lar point changes decisively (see Fig. 0 . The singular point 
is no longer hyperbolic because all linear terms of the planar 
system vanish identically. Instead of that it becomes a gen¬ 
uinely non-linear node with a single characteristic direction 
along the x-axis. In addition there is another singular line at 
x = — | where the derivative of x with respect to In £ tends 
to infinity. 

We must draw our attention again to the upper right 
quadrant where x > — 4 and y > 0. There exist two families 
of integral curves which differ in there asymptotic behavior 
when approaching the singular line at x = — |. Solutions 
that belong to the first family reach the singular line at 
a finite y-value and those from the other family approach 
the singular point at (— |, 0). The latter reside inside the 
shaded area in Fig. 0 These are the only solutions which 
are permitted. 


x 


solutions — 

characteristic direction - 

singular line «***»" 
separatrix — 
non-linear node # 



Figure 5. Phase diagram around the singular point (— §,0) for 
LP viscosity prescription with k = — 1; values on the y-axis are 
scaled by a factor of 10 5 . 


In contrast to the linear analysis discussed above the 
separatrix between these curves is no longer a straight line 
- not even in the immediate vicinity of the singular point. Its 
precise progression is unknown and cannot be deduced an¬ 
alytically, because this would involve solving the non-linear 
problem. Nevertheless one can obtain some important infor¬ 
mation from the non-linear analysifl First of all there exists 
again an unique solution (dashed line in Fig. 0 which ap¬ 
proaches the singular point along the line given by the cubic 
function 

y( x ) = —(* +1) 3 • ( 73 ) 

K + 2 

In the vicinity of the singular point all solutions which pass 
any point above this curve will not proceed to the singular 
point. Instead they will end somewhere on the singular line 
as x — > — | with y > 0. Secondly, all integral curves between 
the unique solution and the characteristic direction y = 0 
approach the singular point along a curve with 

y{x) oc (x+ |) 5 . 

There are infinitely many solutions with this asymptotic be¬ 
havior differing only in their constant of proportionality. 

If we take these results and insert them into the bound¬ 
ary condition Eq. (1681) in the same way as it was done in the 
linear case, one finds that the exponent | is compatible with 
the finite torque boundary condition whereas the power law 
in Eq. (23} would cause the torque to vanish at the inner 
boundary. 

In summary, it can be stated that a feasible solution 
y(x ) of the self-similar disc problem is uniquely determined 
by the parameter k, and the torque applied to the disc at 
the inner boundary. The remaining problem is to specify 
the constant of integration in Eq. (1581) which uniquely de¬ 
termines y as a function of the similarity parameter £ and 
hence D(r, r). This can be achieved in principle by fixing the 


7 We omit the details here and refer to [Frommei] dl92gll for an 
elaborate discussion of the problem. 
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Figure 6. Log-log plot of the two families of integral curves for 
LP viscosity prescription with re — — 1; the dashed line shows the 
unique solution obtained for the no-torque boundary condition. 


DSB/RZ 

LP 

zero torque 1 

3 

finite torque | 

9 

2 


Table 1. Value of the exponent q for different viscosity prescrip¬ 
tions and boundary conditions. 


constant related to the initial condition lEa. 1631) . Unfortu¬ 
nately this would imply that we have to fix the asymptotic 
behaviour as £ —>■ oo in addition to the condition already ap¬ 
plied at the inner boundary as £ —> 0. This would turn the 
initial value problem into a boundary value problem which 
is slightly more difficult to treat numerically. Instead we pro¬ 
ceed in a different way which allows us to apply all auxiliary 
conditions at the inner boundary. 


3.3.2 Solution procedure 


It was already mentioned above that all feasible solutions 
approach the singular point at x = — |. These solutions 
converge towards 

y(t) = yo^ (74) 

as £ —>■ 0 where yo is some not yet determined constant 
which is related to the non-dimensionalized central mass 


M*(t) = lim M = lim r 3 U 2 . 

r—>-0 r—¥ 0 

If we substitute r and Q using the group invariants (Eq. 1531) 
this could be written in terms of £ and y(£) according to 


M*(r) 


r~(i+ 2 ) D 

-2- llm 

re 2 5->0 y * 


r-(S+ 2 ) 

~*vT 


(75) 


where we used Eq. ca to compute the limit. This remark¬ 
able result is not only an analytic formula for the exact 
temporal evolution of the central mass, but allows us also 
to compute the unknown constant yo once we specify the 
central mass M* at some time to. 

In the previous section we showed that in the limit x —► 


— | the function y depends on x according to 

y(x) = y 1 (x+l) q (76) 

where the exponent q takes different values depending on 
viscosity prescription and boundary condition (see Tab.Q]). 
The constant y\ has to be determined from the inner bound¬ 
ary condition. In case of vanishing torque we already derived 
the values for y\ in Eqs. m and (El for all three viscosi¬ 
ties. For the finite torque boundary condition one has to 
compute the limit of Eq. (|671) as r —» 0 which yields 

G*(t) = lim G(r,r) = ~^yf^yf^T~^ + A (77) 

where the constant factor ( depends on viscosity (3 for DSB; 
| for RZ; 12 for LP). We can now compute the value of y i, if 
we specify a distinct value for the torque G* on the inner rim 
at some time to provided that we have already determined 
yo from condition El- 

Summing up the solution procedure there are basically 
six steps: 

(i) select the viscosity parametrization 

(ii) specify /3, re, t 0 , M*(t 0 ), G*(t 0 ) 

(iii) compute yo and j/i 

(iv) choose ^oCl close to the singular point 

(v) compute y(£ 0 ) and x{y(£ o)) 

(vi) integrate the ODE. 

If the torque G*(to) vanishes, the viscous coupling constant 
P affects the solution only indirectly through the time vari¬ 
able t (see Eq. 1501) . The value of £o is somewhat arbitrary 
as long as x(£o) ~ —3/2 holds. It is always possible to scale 
down £o until its impact on the solution is smaller than nu¬ 
merical errors introduced by the integration scheme. 


4 RESULTS AND DISCUSSION 

Before we start the discussion of the numerical results we re¬ 
call some important properties of similarity solutions. The 
requirement that a solution evolves self-similarly implies 
that there exists a relation between its time variable t and 
spatial variable r which becomes manifest in the definition 
of the similarity variable £(r, t) 1 Eg. 1531) . Therefore the func¬ 
tions *(£) and y(£) carry information on both: temporal evo¬ 
lution and spatial dependency. For a fixed instant of time 
£ is simply proportional to the radial coordinate r whereas 
for a constant radius £ is proportional to some power of the 
time variable t. If one uses logarithmic scaling this rela¬ 
tion becomes a simple linear map between the logarithms: 
In Z = lnr + re _1 lnT. Hence we may convert the depen¬ 
dence on In £ into either spatial dependence on In r applying 
a constant shift or temporal evolution by scaling with re 
plus adding a constant shift. Recall that re < 0. Therefore 
increasing £ corresponds to decreasing t and vice versa. 


4.1 Similarity solutions 

In contrast to the rather complex structures shown in the 
phase diagram in Fig.[3]the solutions look quite simple. The 
shape of x(£) is basically a step function with a smooth 
transition between two or three levels depending on the in¬ 
ner boundary condition. If the torque applied to the disc 
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in ? 


Figure 7. Zero torque solutions for different viscosity prescrip¬ 
tions obtained with re = —l,ro = l,M*(ro) = l,ln?o = —10. 
The figure shows the results for both functions a;(?) (bottom left 
to top right, left axis) and y(£) (top left to bottom right, right 
axis). 


at the inner boundary is constant (Fig. [8] upper panel) or 
zero (Fig. Q the two levels are X-^a = —3/2 and *+oo = k. 
In the two cases where the torque at the inner boundary 
changes with time there is an additional intermediate level 
at Xo = —5/4 (Fig. [H] lower panels). Since —x is the log¬ 
arithmic derivative of y with respect to ? (Eq. 15611 y{£) is 
essentially a broken power law. 

Apart from the influence of the inner boundary condi¬ 
tion the solutions show a more or less prominent dependence 
on the viscosity law. The results for DSB and RZ viscosity 
are very similar as one would expect, because the viscos¬ 
ity prescriptions differ by a factor of |a:| which is always of 
order one and becomes at most 3/2 in the Keplerian limit 
where both solutions approach i_ ro = —3/2. In both cases 
the transition between the different constant levels of x is 
quite sharp. In contrast to that the results for x(F) with 
LP viscosity exhibit always a smoother transition. In case 
of decreasing and increasing torque (Fig.[8]lower panels) the 
intermediate level at xo = — 5/4 is hardly visible. 

Once the similarity problem has been solved for the 
non-dimensional functions x(£) and ?/(?) one can extract 
from these the complete information about the self-similar 
evolution using the definition of the group invariants in 
Eqs. (1531) . Thus with Eq. (1541) one obtains as a function of 
r and r and with Eqs. and (EH) one computes the non- 
dimensionalizec(f| expressions for enclosed mass and surface 
density 

M = k~ 2 r~d +2 ') y~ 2 (78) 

£ = ±K~ 2 T-(i +2 ) (2x + 3)Zy~ 2 . (79) 


The local accretion rate can be derived from (1781) by differ¬ 
entiation 


M = d t M = /3d T M = 2/3 k~ 3 +3 ) (x - k) ? 3 y ~ 2 (80) 

and with help of Eqs. © and © one obtains the radial 
velocity 


Up 


M 

2'ktYi 


—2 /3kT 1 t~(» +1 )? 


X — K, 

2x + 3’ 


(81) 


see remark in the first paragraph of Sec. [3] 
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Figure 8. Finite torque solutions for different viscosity prescrip¬ 
tions obtained with to = 1,M*(tq) = l,/3 = 10“ 3 ,G*(ro) = 
— l,ln?o = —25. The figure shows the results for both functions 
ir(?) (bottom left to top right, left axis) and j/(?) (top left to 
bottom right, right axis). 


if k ^ —3/2. The special case where re = —3/2 is discussed 
in more detail below. 

The viscous torque is given by Eq. ©. In addition we 
can derive an ex pression for the vertica lly integr ated dissi¬ 
pation rate (see IKato, Fukue fe Mineshigg |2008i . Chap. 3) 
using Eq. © to substitute the viscosity 

Qvis = uT, ( rdrQ ,) 2 = /3r 2 n 3 -£x 2 f(x). (82) 

We can again utilize the group invariants (1531) , insert E from 
Eq. (l79ll and substitute f(x) with help of Eq. (1571) : 

Qvis = £(-k)~ 5 t~(k + 5 ') xz{x)fy~ 5 . (83) 

Observe that all expressions derived above depend on ?(r, t), 
explicitly as well as implicitly through *(?) and y(£)- Thus in 
addition to the explicit dependence on r there is an implicit 
dependence through the similarity variable ?. 

4.2 Self-similar time evolution 

In figures 0 and © we show the numerical results of the 
k = — 1 solution obtained with zero torque boundary condi¬ 
tion and DSB viscosity prescription. At any fixed time the 
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Figure 9. Self-similar solutions of angular velocity, surface den¬ 
sity and enclosed mass at three different times. All quantities are 
given in non-dimensional units. The figures show results for the 
zero torque boundary condition obtained with DSB viscosity pre¬ 
scription and parameters re = —1 ,/3 = 10" 3 , tq = l,M*(ro) 

■ = - 10 . 


radial dependence of any function is basically a broken power 
law. The kink thereby separates the region in which the cen¬ 
tral object dominates the gravitational potential from that 
in which self-gravity of the disc becomes important. One can 
obtain the asymptotic exponents of these power laws analyt¬ 
ically from the expressions given above using the asymptotic 
expansion of the functions x (£) and y (£). The exponents are 
listed in Tab. [2] 

In addition to the spatial dependence we show the nu¬ 
merical results for three different values of the time vari¬ 
able r in the diagrams. The temporal evolution is given by 
a simple shift in the log-log diagrams. Thereby the radial 
shift is determined by the definition of the similarity vari¬ 
able £ (Eq. 1531) which relates time dependence to spatial 
dependence. Hence the kink where the power law exponent 
changes is shifted by a factor of — ilog(^) on the loga¬ 
rithmic scale if time progresses from n to r 2 . In the limit 
r —> 0 we can compute analytical expressions for the time 
dependence which can again be expressed in terms of power 
laws. The exponents are listed in Tab. [3] In the limit of large 
radii it is not surprising that there is obviously no tempo¬ 
ral evolution because the time scales become larger and the 
solutions do not evolve at all. 




log r 

Figure 10. Self-similar solutions of accretion rate M and radial 
velocity at three different times. All quantities are given in non- 
dimensional units. The figures show results for the zero torque 
boundary condition obtained with DSB viscosity prescription and 
parameters re = —1,/3 = 10“ 3 ,ro = l,M*(ro) = l,ln£o = —10. 

We would like to emphasize that except for surface 
density and radial velocity the power law exponents of the 
asymptotic expansions do not depend on the viscosity pre¬ 
scription. Furthermore angular velocity, enclosed mass and 
accretion rate show always the same asymptotic progres¬ 
sion no matter what kind of viscosity prescription or inner 
boundary condition has been selected. In the Keplerian limit 
enclosed mass and accretion rate become almost constant 
with respect to radial distance. Thus we retain an impor- 
tant feature of the standard steady disc models (see, e.g. 
|Pringlelfl98lh in this limit. However, since we account for 
mass accretion on to the central object, it inevitably be¬ 
comes more massive as time progresses. Therefore the cen¬ 
tral mass as well as the accretion rate generally depend on 
time. There are two exceptions namely the solutions ob¬ 
tained for re = —3/2 and re = —1 which we discuss below. 

In the standard theory of steady state accretion discs 
the viscous dissipation rate Q v i s is proportional to the cen¬ 
tral mass mu ltiplied by th e accretion rate and divided by 
r 3 (see, e.g. IPringlel llQSlI I 3 !. The viscous dissipation rate 
obtained for our solutions with zero torque inner boundary 
condition shows exactly the same asymptotic progress for 
small radii. This becomes clearer if we express Q v is in terms 
of accretion rate and enclosed mass. With help of Eqs. (SU), 
63 and (1801) one can express the product of viscosity and 
surface density in terms of the accretion rate and the simi¬ 
larity solution x((;) and ?/(£): 

U Y, = -AL + 3 )/( j ) 

47 t (x — n)y 


9 There is usually another multiplicative factor due to boundary 
conditions depending on r as well which we ignore here. 
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' ifK = — 3/2 DSB/RZ: exponential decay, LP: no solution 
^ if k = —3/2 zero torque: v r = 0, finite torque: see Eq. d851 
* ifK = —5/4 exponential decay 


Table 2. Asymptotic behaviour of angular velocity, surface den¬ 
sity, enclosed mass, accretion rate, radial velocity, viscous torque 
and dissipation rate for different inner boundary conditions and 
viscosity prescriptions. The table lists the power law exponents 
of the radial dependence for small and large radii. 
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' if k = —3/2 M~k is constant and therefore M vanishes for r 0 
ifK = —3/2 zero torque: v r = 0, finite torque: see Eq. 1 85 1 > 

Table 3. Asymptotic behaviour of angular velocity, surface den¬ 
sity, enclosed mass, accretion rate, radial velocity, viscous torque 
and dissipation rate for different inner boundary conditions and 
viscosity prescriptions. The table lists the power law exponents 
of time evolution in the limit r —> 0. 


If one substitutes v E in Eq. (1821) and uses the balance law 
(f30l) to eliminate fl 2 from the equation, we can express the 
dissipation rate in terms of accretion rate, enclosed mass 
and the similarity solution x(£) and y(£) (recall that we 
use non-dimensional functions and therefore set the gravi¬ 
tational constant G = 1): 


MM xz(x) 
4-7Tr 3 (x — K,)y 


(84) 


Thereby we utilized the viscosity dependent function z(x) 
defined in Eq. ([ST]). One can easily proof in case of the 
zero torque boundary condition that the second fraction in 
Eq. (l84l) becomes 



xz(x) 
(x - n)y 


-3 


in the Keplerian limit for all three viscosity prescriptions. 
Thus we can recover an important relation known from 
steady state accretion disc theory. However, in our case M 
and M are time dependent functions in the limit r —> 0 
(see Tab. [3]) and therefore Q v i s is also time dependent. If we 
examine solutions obtained with the finite torque boundary 
condition we yield slightly steeper radial slopes (see Tab. [5]) 
for the dissipation rate if r —¥ 0. 

In addition we observe a considerable deviation from 
the r~ 3 law depending on the value of k at large radii where 
the disc becomes fully self-gravitating. This is not a surpris¬ 
ing result, because the material in the disc contributes to 
the gravitational potential energy. Hence there exists an ad¬ 
ditional energy source which alters the radial dependence of 
the viscous dissipation rate. 

One easily verifies that the numerical results shown in 
Fig. l9l and [TOl for the k = — 1 solution satisfy exactly the 
predicted asymptotic behaviour listed in Tab. [2] and [3] E.g. 
with those power law exponents we conclude that M oc r°r 1 
for r —> 0. Thus M becomes independent of r for small 
radii resembling the fact that M approaches the value of 
the central mass M* which grows linear with respect to the 
time variable r. Therefore it increases by a factor of 10 3 if 
r increases by the same factor (see Fig. [9] lower panel). 


An exceptional feature of the k = — 1 solution is the 
almost constant accretion rate M (see Fig. 1101 upper panel). 
This is always fulfilled in the limit r —> 0 independent of the 
actual value of k whereas only for k = — 1 one obtains also 
a constant accretion rate for r —» oo (see Tab. [2j. Although 
the accretion rate is almost constant with respect to both 
time and radial distance the solution is not stationary. All 
other quantities show a clear time dependence. 

Another interesting example is that obtained for k = 
—3/2 which has the unique property that the mass of the 
disc is finite and time-independent (see Sec. 13.2.31) . As in 
the other cases there are two different classes of solutioni 33 ! 
depending on the torque applied to the disc at the inner 
boundary. If the torque vanishes one obtains a pure Keple¬ 
rian rotation law with x = —3/2 everywhere at any time. 
Hence the enclosed mass M vanishes except for the central 
mass M*. However, if the torque at the inner boundary is 
finite, there exists another solution with Keplerian rotation 
at small radii as well as larger radii and an intermediate 
flatter rotation law where x approaches —5/4 (see Fig. [TT] 
upper panel). In both cases there is no accretion on to the 
central object. In the former case there is no mass flow at 
all whereas in the latter case the torque applied to the discs 
inner boundary causes radial outflow. The radial velocity 
given by Eq. m with k = —3/2 is positive for all radii 
r > 0: 

Vr = y- = \- t (85) 

The surface density depicted in Fig.[TT]for the solution with 
constant disc mass and decreasing torque shows an expo¬ 
nential decajO at large radii in contrast to all other self¬ 
similar solutions which decline according to a power law 
(see Tab. 0. Hence the disc has a sharp outer rim which 
moves further outwards as time increases. The disc mass is 
redistributed yielding a dispersion of the hole disc. The same 


10 Both solutions are prohibited by the LP viscosity prescription, 
because there is a singular line at x = —3/2 in the phase diagram 
(see Fig. 0 

11 This might be a problem, because it contradicts one of our 
basic assumptions (see Sec. 12.21) . 
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Figure 11. Self-similar time evolution of angular velocity, sur¬ 
face density and enclosed mass for the constant disc mass solu¬ 
tion. All quantities are given in non-dimensional units. The figures 
show results for decreasing torque boundary condition obtained 
with DSB viscosity prescription and parameters k, = — 3/2,/3 = 
10 —3 , r 0 = 1 ,M*(r 0 ) = 1, G* (t 0 ) = -l,ln£ 0 = -20. 

behaviour was also found bv lLin fe Pringlel (ll987l) if they as¬ 
sume in their model, which does not account for self-gravity 
that the disc mass is time-independent. 

4.3 The impact of the power law exponent re 

In Sec. 13.21 we showed that re is the power law exponent 
of the rotation curve at large radii. Since f i is related to 
the enclosed mass (Eq. 130ll whose gradient determines the 
surface density distribution (Eg. 1361) . re also controls the ra¬ 
dial behaviour of these at large radii (see Tab. 0. A steeper 
rotation law causes a steeper radial decline of the surface 
density which in turn leads to a flatter radial increase of the 
enclosed mass. Thus re is just another measure of the matter 
distribution within the disc and hence how self-gravitating 
the disc becomes at large radii. 

On the other hand we found that re has an impact on 
the temporal evolution of the disc at small radii as well. 
From the values of Tab. [3] we can draw some remarkable 
conclusions. First, the central mass always increases with 
time except for the unique solution with re = —3/2 where 
the central mass remains constant. Second, the accretion 



0 2 4 

log r 


Figure 12. Accretion rate M and radial velocity v T as a func¬ 
tion of radial distance r at time t = tq for different values of the 
parameter re. The values on the y-axis are scaled by a factor of 
0~ 1 = 1000. The figures show results for zero torque boundary 
condition obtained with DSB viscosity prescription and parame¬ 
ters 0 = 10 —3 ,ro = l,M*(ro) = l,ln£o = —10. 

rate at small radii decreases with time if re < —1, remains 
constant during the whole evolution only if re = —1 and 
increases with time if re > —1. The slope of the temporal 
decline becomes flatter as re approaches the value —1. Hence, 
discs with higher values of re evolve faster than those for 
which re is close to the Keplerian value of —3/2 or thinking 
in terms of mass accretion rates: Objects embedded in self- 
gravitating discs with flatter rotation laws grow faster than 
those embedded in nearly Keplerian discs. 

If we take a look at the radial dependence of the ac¬ 
cretion rate (Fig. 1121) . we can identify another difference be¬ 
tween solutions obtained with different values of the param¬ 
eter re. We already mentioned before that only for re = — 1 
the accretion rate becomes independent of radial distance 
in the limit r —> oo. If re is above that value, the accretion 
rate increases as r increases, it tends to zero if re < —1 and 
it remains positive for all radii if re ^ —5/4. However, if 
re < —5/4, the accretion rate falls below zero in the transi¬ 
tional region between the inner Keplerian disc and the outer 
self-gravitating part of the disc. The same behaviour can be 
seen more clearly if we take a look at the radial velocity 
depicted in the lower panel of Fig. 1121 For small radii the 
radial velocity is always negative leading to the positive ac- 
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cretion rate in the region where the disc is nearly Keplerian, 
but as r increases the radial velocity becomes positive for 
the solutions obtained with k = —1.4 and k = —1.3. These 
accretion discs may loose a considerable amount of mass by 
redistributing it to the outer disc rather than accreting it on 
to the central object. This supports the proposition stated 
in the previous paragraph. 


4.4 Verification of the model assumptions 

The derivation of the disc evolution equation in Sec. m 
relies on two essential assumptions, namely the thin disc 
approximation which is related to the supersonic rotation 
requirement (see Eq. 12.41) and the slow accretion limit. The 
verification of the former would require knowledge about 
the speed of sound and therefore the temperature in the 
midplane of the disc (see Eq. m- Hence in order to check for 
supersonic rotation, one generally has to solve the vertical 
structure problem which is beyond the scope of this workF 2 ! 

Nevertheless we can compute the ratio v r /v v and figure 
out if the radial velocity is always much smaller than the 
azimuthal velocity. This is at least a necessary condition for 
the slow accretion limit. From Eqs. (El and m one easily 
computes 


«r = ^ (a- - k) y 
V(p x T ~2 


( 86 ) 


Hence v r /v v scales with the non-dimensional viscous cou¬ 
pling constant f) and depends on the non-dimensional func¬ 
tions x(£,) and y(£). We can use the numerically obtained 
similarity solutions to compute the inverse £(a:) and insert 
the result in y(£) to express the right hand side of Eq. (1861) 
in terms of x. This allows us to express the velocity ratio as 
a function of x alone. 

The results for various numerical solutions obtained 
with different values of the parameter k and with different 
viscosity prescriptions and boundary conditions are depicted 
in Fig. 1131 One can clearly see that in all cases v r /v v /3 _1 is 
of the order of one and therefore v r /v v must be of the order 
of p which is usually set to values smaller than 10 -2 (see 
Sec. 12.611 . Therefore we conclude that as long as P <S 1, the 
radial velocity v r is indeed much smaller than the azimuthal 
velocity v v and furthermore if c s , c is only slightly larger than 
v r the azimuthal velocity must be highly supersonic. 


4.5 Application to AGN discs 

In this section we exemplify how to introduce physical 
dimensions and apply them to the non-dimensional solu¬ 
tions shown in the previous sections. As was already men¬ 
tioned in the introductory paragraph of Sec. El the non- 
dimensionalization of the disc equations can be achieved by 
specifying two independent basic scales, e.g. , length and 
mass. Then a third independent scale, e.g. , time, is given 
by Eq. ED- One can specify any two of the three scales 
for length, mass and time and compute the remaining scale. 
Once the basic scales have been determined one can derive 


12 We explain in Section 14. 5 1 how to compute a rough estimate of 
the central temperature for a typical AGN accretion disc. 
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Figure 13. Ratio of radial and azimuthal velocity v r /v<p as a 
function of local power law exponent x. The top panel shows 
results obtained with DSB viscosity and zero torque boundary 
condition for different values of the parameter k. The lower panel 
shows solutions for different viscosity prescriptions and boundary 
conditions for k = — 1. The values on the y-axis are scaled by 
a factor of f3~ 1 = 1000. The other simulation parameters are 
ro = 1,M*(to) = l,ln£o = —10 for zero torque solutions and 
ln£o = —25 for finite torque solutions with G*(ro) = — 1. 


the scales for angular velocity and linear velocity, surface 
density and viscous dissipation rate: 

Q = r -1 , v = ff _1 , E = Mr -2 , Q = Mr ~ 3 

and since r = fit one can compute the scale of the real time 
variable t and the mass accretion rate 

t = j8 _1 f, M = Mr 1 . 

In Tab. E] we list the scaling constants of a self-gravitating 
accretion disc model applicable to AGN. These scales can be 
used as units for the non-dimensional model shown in Fig. El 
e.g. the numbers on the x-axis denote the common logarithm 
of radial distance to the central black hole in astronomical 
units (AU). To convert the radial scale from AU to parsec 
(pc) one has to subtract 5.3. Hence the outer boundary at 
10 8 AU corresponds to 10 2 ' 7 pc ~ 500 pc or 1.5 ■ 10 21 cm. If 
we set p = 10 -3 (for a discussion of reasonable values for /3 
see Sec. 12. Gl and the references given there) we can determine 
the constant accretion rate of the central black hole from the 
upper panel of Fig. [101 

M = pM = 10 -3 • 200 M 0 yr _1 s = 0.2 M 0 yr _1 . 
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r 

f 

t Q 

V 

1 AU 

5 ■ 10 -3 yr 

5yr 200 yr —1 

948 kms -1 

M 

£ 

M 

Q 

10 3 Mq 

9 ■ 10 9 g cm -2 

200 Mq yr —1 5 

■ 10 17 Wm- 2 


Table 4. Scales for the AGN example with 0 = 10 —3 . Most of 
the values have been rounded to one significant digit to fit into 
the table. For typical paramete r s of an evolved AGN di sc see 
Colli n-Souffrin fe D umont) (i 199f| j: [bin &; Papaloizo~ul Jl996T) ■ 

This is a quite reasonable value for the late phases of 
the accretion process. However, if the mass of the central 
black hole is below 10 7 Mq this is well above the maxi¬ 
mum accretion rate permitted by the Eddington limit which 
is - assuming an accretion efficiency of 10% - of the or¬ 
der of M max ~ 2 ■ 10 _8 (M/Mq) Mq yr _1 . Nevertheless, if 
one assumes that the black hole accretes at the Edding¬ 
ton limit during its early evolution switching to the con¬ 
stant sub-Eddington accretion rate after the Salpeter time 
scale JSalpeterlll964T ). it is still possib le to accumulate up 
to 10 9 Mq within a few billion years (iDuschl fc Strittmatted 

B). 

If one assumes that energy transfer inside the disc is 
mainly due to radiation in the vertical direction one can 
estimate the effective temperature of the disc by equating 
the energy generation due to viscous dissipation and energy 
loss caused by radiative cooling: 

Qvis = Qcool = 2(7 Teff (87) 

where a is the Stefan-Boltzmann constant. The factor 2 is 
necessary to account for both radiating surfaces of the disc. 
Since Q v i s is completely determined by the self-similar solu¬ 
tion (see Eg. 1831) . we can solve the equation above for T e g. 

The result for the AGN example is shown in Fig. [TT] (up¬ 
per panel). In principle we can compute the temperature dis¬ 
tribution for arbitrary small radii, but below a certain mini¬ 
mal radius it becomes meaningless. Therefore we decided to 
truncate the curves at the last stable circular orbit around a 
Schwarzschild black hole l ocated at three Schwarzschild radii 
dNovikov fe Thornelll973h . In view of the previously shown 
solutions, it is not astonishing that the radial temperature 
profile is given by a broken power law. The exponents are de¬ 
termined by the exponents of Q v is (see Tab. 0 multiplied by 
a factor of 1/4 which becomes —3/4 at small radii in case of 
the zero torque boundary condition and (5k + 3)/4 = —1/2 
(with re = —1) in the self-gravitating outer regions of the 
disc. The slope of the radial temperature profile, say p, is 
related to the slope of the s pectral energy dist ribution ac¬ 
cording to n = 3 — 2/p Csee lLvnden-BellM1969l ). Again our 
model recovers the spectral index n = 1/3 for Keplerian 
discs, but in the self-gravitating region this becomes n = — 1 
leading to an enhanced radiative flux at higher wave lengths. 

So far we have avoided the intricate computation of the 
vertical disc structure. This would in general require the 
solution of the energy equation in conjunction with the ra¬ 
diative transfer equation which is beyond the scope of this 
work. However, if we assume an isothermal vertical structure 
- which is a very gross simplification - the midplane tem¬ 
perature equals the effective temperature. This allows us to 
estimate the speed of sound and the azimuthal Mach number 
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Figure 14. Self-similar time evolution of a massive AGN disc. 
The diagrams show effective temperature T e g, Mach number 
Vtp/cs : c, scale height h and Toomre parameter Q as a function of 
radius at three different times. The results were obtained with the 
same parameters as those depicted in Figs. iQlandlTOl namely zero 
torque boundary condition with DSB viscosity prescription and 
parameters re = —1,0 = 10“ 3 ,ro = l,M*(ro) = l,ln£o = —10. 
The solutions are truncated at 3 Schwarzschild radii (see text). 


v v /cs, c in the equatorial plane depicted in the second panel 
of Fig. [14] The results show that the rotational velocity is 
indeed supersonic which is required by the model assump¬ 
tions (see Sec. 12. 41) . In the early phases of the AGN evolution 
the flow is only slightly supersonic with a minimum Mach 
number of the order of 10 around the region between the Ke¬ 
plerian and the self-gravitating disc. The minimum remains 
in this region during the evolution of the AGN, but its value 
increases up to roughly 10 3 . 

With the estimate of the azimuthal Mach number one 
can verify another important requirement for the model, 
namely the thin disc assumption h/r -C 1. This ratio of 
scale height and radius is shown in the third diagram of 
Fig. 1141 Since this ratio is related to the azimuthal Mach 
number according to Eq. the assumption is again bet¬ 
ter for the evolved AGN disc. However, it remains always 
smaller than one but reaches values of the order of ICG 2 in 
the early phases of the evolution. 

Although our cooling model is very simple, our results 
on effective temperature, Mach number and aspect ratio are 
in good agreement with those given in the literature (see e.g. 

ICollin-Souffrin fe Dumon3ll990l : lLin &; Papaloizoulll996l ). at 

least in case of an evolved AGN disc. Besides that one should 
be aware of the fact that any results on the early evolu- 
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tion of an AGN disc are quite speculative, because obser¬ 
vational constraints on accretion discs around intermediate 
mass black holes with masses in the range of 10 3 to 10 s Mg 
are very limited. Right n ow t here seems to be one promising 
candidate dDavis et alj|201lh . 

Another important quantity when talking about self- 
gravitating accretion discs is the Toomre stability parame¬ 
ter Q defined in Eq. (14411 . If its value is below one, the disc 
becomes gravitationally unstable. Keplerian discs are always 
stable. This can be seen most easily with help of Eq. (1451) . 
If the logarithmic gradient x of the rotational velocity ap¬ 
proaches the Keplerian value Q becomes infinitely large as 
long as the inverse of the azimuthal Mach number c %tC /v v 
remains finite. On the other hand, if x deviates from its Ke¬ 
plerian value of —3/2 the disc may become unstable if the 
rotational velocity is highly supersonic. In Fig. [44] (lower 
panel) the Toomre parameter for the AGN disc is shown as 
a function of radius at three different times. The outer re¬ 
gions of the disc are always unstable, because the slope of 
the rotational velocity deviates considerably from its Kep¬ 
lerian value and the Mach number becomes rather high. In 
the early phases of the evolution the transition from sta¬ 
ble to unstable occurs roughly in the intermediate region 
where the disc becomes fully self-gravitating. This location 
is shifted towards lower radii as time progresses, because the 
Mach numbers in the evolved discs are higher even in the 
region where the discs are almost Keplerian. 

One might argue that the existence of these solutions is 
questionable, because the gravitational instabilities would 
lead to fragmentation of the whole disc. Since the time 
scales for these local gravitational collapses are rather short 
compared to the viscous time scale, the disc would disinte¬ 
grate into small clumps forming stars before a considerable 
amo unt of gas has been accreted on to the c entral black 
hole dLin fc Pringlell 19871 : IShlosman et aLlll990l ’). Since these 
processes are v ery sensitive to the temperature of the disc 
dGammiel 1200il l one should be careful with this reasoning, 
because we applied a far too simple cooling model to de¬ 
rive the temperature distribution. As was already mentioned 
above, the deduction of a reliable disc temperature is quite 
difficult and beyond the scope of this paper. Therefore we 
leave it for future work. 

4.6 Comparison with other works 

The literature on one-dimensional modelling of self- 
gravitating accretions discs using the local viscous trans¬ 
port approximation can be subdivided into two different 
categories: (i) Keplerian models in which the rotation law 
is not affected by self-gravity and (ii) fully self-gravitating 
disc models in which one considers the discs potential in the 
radial balance law. In both cases there exist stationary as 
well as non-stationary models. We briefly discuss the for¬ 
mer here, because they are not subject of the present work, 
and compare our results in more detail with other time- 
depende nt m odels. 

IPaczvnsldl (119T8I 4 proposed the first self-gravitating 
AGN model considering an evolved AGN with supermas- 
sive black hole mass of 10 10 M@ surrounded by a rather 
compact Keplerian disc (Mdisc = 10 9 M@). The surface den¬ 
sity in this model shows a steep decline (oc r~ 3 ), which we 
never yield in our models due to the constraints on re (see 


Tab. 

2|. Other authors ( 

Mineshise & Umemura 19961; Bert in 

1997 

Bcrtin & Lodato 

1999i: Duschl et al.1 20004 find fullv 


self-gravitating stationary solutions with different viscosity 
prescriptions including a-viscosity (with self-regulation) and 
/3-viscosity. All these models show a remarkable common fea¬ 
ture, namely that the rotation curves become flat and that 
E oc r _1 in the limit r — » oo. We find the same asymptotic 
behaviour and an almost constant accretion rate in our self¬ 
similar quasi-stationary model which we obtain for re = — 1 
(see Sec. 14.21) . 

iLin fe Pringle! (1 19871 4 developed a time-dependent self¬ 
similar model for self-regulated Keplerian discs accounting 
for self-gravity only in their viscosity model (see Sec. 12.61) . 
They found an analytic solution for the zero torque bound¬ 
ary condition if total disc angular momentum is conserved. 
Within the context of our self-similar model we already 
showed that the case of constant disc angular momentum 
is not permitted by the basic constraints on the slope of the 
rotation curve, because it would require that re = —5/3 (see 
Sec. 13. 2. 3D . The analysis of the phase plane for our model 
reveals that in this particular case the local exponent of the 
rotation law x becomes smaller than the Keplerian value of 
—3/2 at some finite radius for all solutions which are Keple¬ 
rian near the origin. Thus beyond this radius the slope of the 
rotation law is steeper than permitted by physical reasons 
and as a consequence E becomes negat i ve (se e Eq. l37D . Inter¬ 
estingly the solution of lLin fe Pringle! d 19871 4 yields complex 
values for E beyond a certain radius. Hence this solution 
also breaks down at the outer rim of the disc. Apart from 
these general problems our results exhibit some features of 
their solution near the origin, namely that E oc r _3//2 for all 
solutions obtained for LP-viscosity with no central coupling 
(see Tab. [2j) and that M oc t 6 ^ 5 for re = —5/3 (see Tab. [3]). 

iRice fc Armitagel (j2009l 4 extended this model consid¬ 
ering a cooling mechanism which allows to determine the 
a parameter o f the effective viscosity self-consistently (see 
[Gammiell200ll 4. Although the authors claim to solve a fully 
self-gravitating disc model, we would prefer the term Kep¬ 
lerian model here, because they do neither account for the 
disc material in the radial balance law (Eq. I28D nor do they 
consider the impact of self-gravity on the scale height (see 
Sec. 12. H and Eq. l38D . Nevertheless, we can compare their nu¬ 
merical solutions for the time evolution of the surface den¬ 
sity with our self-similar model. A remarkable common fea¬ 
ture is again that E oc r _3,/2 for small radii corresponding 
to our model with LP-viscosity and zero torque boundary 
condition. At about 2AU the slope of the surface density 
becomes steeper (E oc r~ 2 ) which is the asymptotic limit 
as r —> oo for almost massless discs (re m —3/2) in our self¬ 
similar model (see Tab. H. Then at about 20 AU the decline 
in surface density steepens again. According to the authors 
this is due to their cooling model which causes a sudden 
drop in temperature. Since our model does not account for 
any cooling mechanism one cannot expect to reproduce this 
feature. 

Based on their previous work ILin fe Pringlel rtl990l 4 
extended their model accounting for self-gravity (in the 
monopole approximation), radiative cooling and energy 
transport. Furthermore they modified the viscosity prescrip¬ 
tion adding th e standard Sh ak ura- Sunyaev viscosity to the 
original model iLin fc Pringld[l987l 4 to avoid the problem of 
a vanishing viscosity coefficient in the inner Keplerian parts 
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LP-model 

a 

K 

h 

-(3/ 

ft 2) 

At 

6 

5 

11 

10 

0.73 

8 ^ 
11 "" 

: 0.73 

B1 

5 

4 

9 

8 

0.67 

2 ^ 
3 ~ 

0.68 

A4 

13 

10 

23 

20 

0.61 

14 ^ 
23 "" 

: 0.61 

B4 

4 

3 

7 

6 

0.58 

4 ^ 
7 ~ 

0.57 


Table 5. Comparison of radial and temporal power law exponents 
of surface density S oc r a and central mass M* oc I J 1 ' derived from 
the self-similar model with values obtained from iLin fe Pringle! 
dl99(lh (columns 2 and 4). 


of the disc. They performed numerical simulations for differ¬ 
ent initial conditions. The radial profiles of surface density 
evolve self-similar in all their simulations with power law ex¬ 
ponents between —6/5 and —4/3 (depending on the model) 
in the fully self-gravitating region. It flattens close to the 
inner boundary of the disc where the Shakura-Sunyaev vis¬ 
cosity dominates. This is consistent with our observation 
that E oc r -1 / 2 for r —» 0 with DSB /?-viscosity and zero 
torque boundary condition (see Tab. 0. We already showed 
in Sec. l2.6l that /I-viscosity recovers a-viscosity in this limit. 

Furthermore we find a remarkable match comparing 
the time evolution of the central mass with our self-similar 
model. If one reads off from their diagrams the power law 
exponents of the radial surface density distribution (at a 
specific time!), say a, one can compute the associated expo¬ 
nents of the rotation law according toK= (a — 1)/2 (see last 
column of Tab. El). With these values for re one can make a 
prediction for the time evolution of the central mass using 
the asymptotic exponent given in the third row of Tab. E] 
From their model Al, for instance, we get a = —6/5 and 
therefore re = —11/10. Thus we predict that M* oc t 8 ^ 11 . 
This is almost exactly the value one can read off their di¬ 
agram (ss 0.73) for the period in which E evolves quasi¬ 
stationary with a constant radial slope (between t = 2ts 
and 3 iff). The exponent for four different numerical mod¬ 
els and our predictions are shown in Tab. [5] This strongly 
supports our observation that the self-similar time evolu¬ 
tion of self-gravitating discs is mainly controlled by a single 
dimensionless parameter, namely re. 

Fully self-gravitati ng self-similar accretion disc m odels 
were first developed in IMineshig e fe Umemural d 1997h and 
iMineshige et al.l d 1997h . In the first paper the authors as¬ 
sume that the disc is isothermal with respect to r whereas 
in the second paper they extend the model using a poly¬ 
tropic relation. Both papers account for self-gravity in the 
monopole approximation and simplify the equations using 
the slow accretion limit. The main difference compared 
to our model is the effective viscosity which depends lin¬ 
early on radius or - in the second paper - according to 
a power law with a constant exponent. Thus the viscos¬ 
ity does not change as the disc evolves in time. Another 
important difference of these models is the similarity trans- 
formation applied to the syste m of differential equations. In 
IMineshige fc Umemural (Il997l l the similarity variable is pro¬ 
portional to r/t which corresponds to our quasi-stationary 
model obtained with re = — 1. At large radii the asymptotic 
behaviour of these similarity solutions is in compliance with 
our findings (E oc r -1 , 9 oc r -1 , M = const, v r = const). 
Whereas for small radii the model fails, because the rotation 


law does not become Keplerian, instead 9 oc r -4//3 . Further¬ 
more the accretion rate is proportional to r 1,73 in the vicinity 
of the central object. It therefore vanishes in the limit r —> 0. 
The results for the polytropic model do not differ substan¬ 
tially from the isothermal model and suffer from the same 
limitations. 

On the basis of the isothermal model presented in 
IMineshige fc Umemural (Il997f) iTsuribel dl999l ~) developed an 
isothermal self-similar a-disc solution depending on a vari¬ 
able Toomre parameter. He utilized the viscosity prescrip¬ 
tion in Eq. m together with the definition of the Toomre 
parameter (Eg. 1441) . This allows him to express the effective 
viscosity in terms of 9 and E in the same way as we do to 
derive the LP-viscosity model (Eg. 1461) . The only difference 
between both viscosity prescriptions is that we additionally 
assume that the disc is in a marginally stable state, i.e. 
<5 = 1. Thus his similarity solutions depend on the param¬ 
eter Q. Apart from this, his model is identical to our quasi- 
stationary model for LP-viscosity and zero torque boundary 
condition with re = — 1 . The comparison with our results 
reveals that all observables show exactly the same asymp¬ 
totic behaviour (see Tab. Eland [3}. Hence, this solution is a 

special case of our m ore general res ul ts. _ 

_ More recently lAbbassi et all (I2006T ): lAbbassi et all 

(|2013jh proposed self-similar disc models using the 13- 
viscosity prescription. As was already pointed out in the 
introduction of the present work, these models are based 
on contradictory assumptions. Thus consequently the so¬ 
lutions presented in these works violate a fundamental re¬ 
quirement for the applicability of the slow accretion limit: 
hig h ly supersonic azi muthal velocity v v c s , c (see Fig. 2 of 
lAbbassi et al.l ll2006lf '). Hence, we do not comment further 
on their results. 


5 SUMMARY 

In this work we have developed a rather comprehensive dy¬ 
namical model for geometrically thin accretion discs. It ac¬ 
counts for time-dependent angular momentum redistribu¬ 
tion and mass accretion on to the central object as well 
as the impact of self-gravity on the accretion process. Sev¬ 
eral reasonable assumptions have been made to reduce the 
complexity of the problem keeping the model as simple as 
possible including the thin disc approximation, the slow ac¬ 
cretion limit and the monopole approximation. On the basis 
of these assumptions we have deduced a new PDE describing 
the time evolution of the angular velocity profile in viscous, 
self-gravitating accretion discs. This can be seen as the ma¬ 
jor difference between our approach and the standard theory 
where disc dynamics is based on the evolution of surface den¬ 
sity. However, since angular velocity is related to the mass 
distribution we can derive the surface density from our so¬ 
lutions which makes our model coequal, albeit more simple 
in view of the underlying mathematical problem. 

We have discussed three different well-established vis¬ 
cosity models taken from the literature, all of them have 
already been applied to self-gravitating accretion discs. It is 
worthwhile to mention here that although one has to specify 
a viscosity prescription in order to solve the disc evolution 
equation, our model is far more general in the first place, 
because it does not rely on a specific viscosity model. We 
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have furthermore shown that in case of the three viscosity 
prescriptions the disc evolution equation admits similarity 
solutions. The associated similarity transformation depends 
on a single non-dimensional parameter which is related to 
the radial slope of the angular velocity at large radii. We 
have shown that this parameter has a major impact on 
the whole disc evolution and that flatter rotation laws yield 
higher accretion rates. Thus we conclude that the more self- 
gravitating the accretion disc, the faster the growth of the 
central object. 

Finally, we have applied our model to an AGN accre¬ 
tion disc. We have discussed the formation of supermassive 
black holes in this context and found that our model of self- 
gravitating accretion discs may explain the rapid growth of 
SMBHs in the early universe. 
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Figure Al. Dimensionless weight functions 
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APPENDIX A: THE FUNCTIONS T< AND JF> 

The definite integral in Eq. (1321) exists for any combination 
of the parameters r ^ 0 and s is 0 where at least one of them 
is different from zero. An analytical expression in terms of 
hypergeometric functions is given in IGradshtevn fe Rvzhikl 
ll2000l) (Eq. 6.574). It may be rewritten with help of the 
complete elliptic integrals of the first and second kind K(k) 
and E(k) according to 

E<(k) = ^-E(k) if = 1 

E >{ k) = \m^l if 

where k is the elliptic modulus. One easily proves using the 
asymptotic expansions of E(k) and K(k ) that 

T<(k) = 1 - (|) 2 + O ( k A ) and F>{k) = § +O ( k 3 ) . 

In order to compute the integral expressions in Eq. (1331) one 
has to evaluate the weight functions (1 — J-<)/k and 7 r >/fc 
plotted in Fig. IA1I 
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